options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

regression

normal linear model

ex5-1.stan

normal regression

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> s;
}
model{
  vector[N] m=X*b;
  y~normal(m,s);
}
generated quantities{
  vector[N] y1;
  vector[N] m1=X*b;
  for(i in 1:N){
    y1[i]=normal_rng(m1[i],s);
  }
}
n=30
mdl=cmdstan_model('./ex5-1.stan')

single regression

tb=tibble(x=runif(n,0,9),
          y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)

qplot(data=tb,x,y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -115.099 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20      -13.7106   0.000154046   0.000912282      0.9892      0.9892       24    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
mle
##  variable estimate
##    lp__     -13.71
##    b[1]      -0.05
##    b[2]       0.97
##    s          0.96
##    y1[1]      8.70
##    y1[2]      9.13
##    y1[3]      8.41
##    y1[4]      6.23
##    y1[5]      8.22
##    y1[6]      5.17
##    y1[7]      3.04
##    y1[8]      2.88
##    y1[9]      6.78
##    y1[10]     6.76
##    y1[11]     1.72
##    y1[12]     7.52
##    y1[13]     6.35
##    y1[14]     5.34
##    y1[15]     2.95
##    y1[16]     3.87
##    y1[17]     0.13
##    y1[18]     2.68
##    y1[19]     1.97
##    y1[20]     3.53
##    y1[21]     3.88
##    y1[22]     4.25
##    y1[23]     1.04
##    y1[24]    -0.71
##    y1[25]     0.00
##    y1[26]     0.41
##    y1[27]     2.23
##    y1[28]     5.19
##    y1[29]    -0.12
##    y1[30]     2.67
##    m1[1]      8.39
##    m1[2]      8.69
##    m1[3]      6.75
##    m1[4]      7.05
##    m1[5]      7.88
##    m1[6]      5.92
##    m1[7]      3.61
##    m1[8]      2.98
##    m1[9]      5.59
##    m1[10]     8.23
##    m1[11]     1.99
##    m1[12]     4.84
##    m1[13]     5.32
##    m1[14]     5.93
##    m1[15]     2.87
##    m1[16]     3.38
##    m1[17]    -0.03
##    m1[18]     2.60
##    m1[19]     3.58
##    m1[20]     4.37
##    m1[21]     3.77
##    m1[22]     5.93
##    m1[23]     1.89
##    m1[24]     0.17
##    m1[25]     0.94
##    m1[26]     0.57
##    m1[27]     2.28
##    m1[28]     4.86
##    m1[29]     1.39
##    m1[30]     2.46
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.35 -14.99 1.36 1.09 -18.00 -13.90 1.00      697      829
##    b[1]    -0.05  -0.06 0.39 0.37  -0.69   0.59 1.02      315      494
##    b[2]     0.97   0.97 0.08 0.07   0.84   1.10 1.02      392      658
##    s        1.04   1.02 0.14 0.13   0.83   1.29 1.00     1013      769
##    y1[1]    8.37   8.35 1.13 1.09   6.56  10.20 1.00     1836     1661
##    y1[2]    8.69   8.70 1.14 1.10   6.80  10.52 1.00     1979     2013
##    y1[3]    6.77   6.76 1.07 1.04   5.04   8.49 1.00     1799     1937
##    y1[4]    7.03   7.04 1.08 1.03   5.26   8.81 1.00     2026     1971
##    y1[5]    7.88   7.90 1.12 1.07   6.06   9.68 1.00     2066     1807
##    y1[6]    5.89   5.91 1.07 1.08   4.08   7.58 1.00     1918     1893
##    y1[7]    3.65   3.66 1.08 1.08   1.90   5.38 1.00     2001     1820
##    y1[8]    3.00   3.00 1.07 1.03   1.23   4.78 1.00     1686     1721
##    y1[9]    5.58   5.59 1.09 1.08   3.75   7.36 1.00     1958     1905
##    y1[10]   8.24   8.23 1.16 1.14   6.32  10.16 1.00     2056     1771
##    y1[11]   2.01   1.97 1.10 1.07   0.23   3.88 1.00     1643     1711
##    y1[12]   4.85   4.86 1.06 1.01   3.14   6.66 1.00     1848     1908
##    y1[13]   5.34   5.35 1.07 1.01   3.56   7.08 1.00     1846     1990
##    y1[14]   5.92   5.92 1.06 1.00   4.13   7.66 1.00     1881     1855
##    y1[15]   2.86   2.86 1.06 1.03   1.17   4.59 1.00     1743     1931
##    y1[16]   3.38   3.40 1.10 1.07   1.56   5.19 1.00     2042     1562
##    y1[17]  -0.05  -0.07 1.15 1.13  -1.85   1.87 1.00     1320     1425
##    y1[18]   2.60   2.59 1.07 1.08   0.88   4.41 1.00     1691     1743
##    y1[19]   3.62   3.59 1.07 1.01   1.87   5.34 1.00     1813     1757
##    y1[20]   4.40   4.40 1.06 1.03   2.62   6.13 1.00     2061     1969
##    y1[21]   3.73   3.73 1.05 1.04   2.03   5.48 1.00     2014     1930
##    y1[22]   5.96   5.98 1.06 1.04   4.22   7.67 1.00     1988     1998
##    y1[23]   1.89   1.86 1.08 1.10   0.08   3.66 1.00     1871     1974
##    y1[24]   0.19   0.18 1.09 1.10  -1.51   1.94 1.00     1443     1837
##    y1[25]   0.91   0.90 1.11 1.10  -0.85   2.74 1.00     1155     1573
##    y1[26]   0.60   0.57 1.09 1.04  -1.25   2.37 1.00     1590     1690
##    y1[27]   2.28   2.29 1.06 1.04   0.57   3.98 1.00     1955     1807
##    y1[28]   4.84   4.84 1.06 1.00   3.17   6.61 1.00     1913     1853
##    y1[29]   1.45   1.45 1.10 1.07  -0.40   3.26 1.00     1650     1646
##    y1[30]   2.45   2.47 1.06 1.04   0.69   4.15 1.00     1899     1793
##    m1[1]    8.39   8.39 0.39 0.37   7.73   9.05 1.01     1003     1238
##    m1[2]    8.69   8.69 0.42 0.38   7.99   9.39 1.01      945     1256
##    m1[3]    6.75   6.74 0.29 0.28   6.28   7.23 1.00     1663     1544
##    m1[4]    7.05   7.05 0.31 0.29   6.54   7.56 1.00     1469     1479
##    m1[5]    7.88   7.88 0.36 0.34   7.29   8.49 1.00     1128     1299
##    m1[6]    5.92   5.92 0.25 0.23   5.52   6.33 1.00     2172     1521
##    m1[7]    3.61   3.61 0.21 0.20   3.27   3.96 1.01      737     1038
##    m1[8]    2.98   2.99 0.22 0.21   2.62   3.35 1.01      501      836
##    m1[9]    5.59   5.59 0.23 0.22   5.22   5.98 1.00     2216     1467
##    m1[10]   8.23   8.23 0.38 0.36   7.59   8.87 1.00     1038     1290
##    m1[11]   1.99   2.00 0.26 0.24   1.55   2.42 1.02      367      647
##    m1[12]   4.84   4.84 0.21 0.20   4.51   5.18 1.00     2037     1613
##    m1[13]   5.32   5.32 0.22 0.21   4.97   5.69 1.00     2202     1592
##    m1[14]   5.93   5.93 0.25 0.23   5.53   6.34 1.00     2171     1521
##    m1[15]   2.87   2.88 0.22 0.21   2.50   3.24 1.01      476      781
##    m1[16]   3.38   3.38 0.21 0.21   3.03   3.73 1.01      625      934
##    m1[17]  -0.03  -0.04 0.38 0.37  -0.67   0.61 1.02      315      489
##    m1[18]   2.60   2.60 0.23 0.22   2.21   2.99 1.01      430      715
##    m1[19]   3.58   3.58 0.21 0.20   3.24   3.92 1.01      718     1027
##    m1[20]   4.37   4.37 0.20 0.20   4.05   4.70 1.00     1491     1522
##    m1[21]   3.76   3.77 0.20 0.20   3.43   4.10 1.00      832     1110
##    m1[22]   5.93   5.93 0.25 0.23   5.52   6.34 1.00     2171     1521
##    m1[23]   1.89   1.89 0.27 0.25   1.44   2.33 1.02      361      638
##    m1[24]   0.17   0.17 0.37 0.35  -0.45   0.78 1.02      316      479
##    m1[25]   0.94   0.94 0.32 0.31   0.40   1.48 1.02      326      549
##    m1[26]   0.57   0.57 0.35 0.33  -0.01   1.14 1.02      320      534
##    m1[27]   2.28   2.28 0.25 0.23   1.86   2.68 1.02      391      670
##    m1[28]   4.86   4.86 0.21 0.21   4.53   5.20 1.00     2046     1612
##    m1[29]   1.39   1.39 0.29 0.28   0.89   1.88 1.02      338      579
##    m1[30]   2.46   2.46 0.24 0.23   2.06   2.85 1.01      411      715

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)

par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

multiple regression

tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
          y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -36.1196 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -13.6854   0.000193542    0.00129532           1           1       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -13.69
##    b[1]       1.40
##    b[2]       0.84
##    b[3]      -1.12
##    s          0.96
##    y1[1]     -1.28
##    y1[2]      2.02
##    y1[3]     -4.52
##    y1[4]      5.49
##    y1[5]     -3.82
##    y1[6]     -2.10
##    y1[7]     -1.27
##    y1[8]      5.65
##    y1[9]     -0.65
##    y1[10]     2.97
##    y1[11]     1.38
##    y1[12]     1.75
##    y1[13]     8.74
##    y1[14]     2.60
##    y1[15]    -2.98
##    y1[16]     3.43
##    y1[17]     0.47
##    y1[18]     6.46
##    y1[19]    -0.90
##    y1[20]     1.03
##    y1[21]     3.67
##    y1[22]     1.96
##    y1[23]    -6.01
##    y1[24]     0.90
##    y1[25]    -0.11
##    y1[26]     3.25
##    y1[27]    -2.32
##    y1[28]    -1.46
##    y1[29]     4.81
##    y1[30]     5.89
##    m1[1]     -1.60
##    m1[2]      0.46
##    m1[3]     -2.91
##    m1[4]      5.70
##    m1[5]     -5.14
##    m1[6]     -1.86
##    m1[7]     -1.59
##    m1[8]      6.31
##    m1[9]     -0.22
##    m1[10]     3.90
##    m1[11]     3.28
##    m1[12]     2.74
##    m1[13]     7.93
##    m1[14]     2.66
##    m1[15]    -2.73
##    m1[16]     2.06
##    m1[17]     1.27
##    m1[18]     6.81
##    m1[19]    -0.90
##    m1[20]     2.91
##    m1[21]     2.36
##    m1[22]     2.05
##    m1[23]    -5.62
##    m1[24]    -0.12
##    m1[25]    -1.31
##    m1[26]     1.37
##    m1[27]    -2.45
##    m1[28]    -1.83
##    m1[29]     6.78
##    m1[30]     4.90
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.90 -15.58 1.49 1.36 -18.91 -14.13 1.01      513      399
##    b[1]     1.41   1.44 0.75 0.77   0.22   2.58 1.01      255      262
##    b[2]     0.84   0.83 0.09 0.09   0.69   0.99 1.02      306      373
##    b[3]    -1.12  -1.12 0.09 0.08  -1.27  -0.98 1.00      510      590
##    s        1.06   1.05 0.15 0.14   0.85   1.34 1.00     1105      941
##    y1[1]   -1.60  -1.60 1.13 1.13  -3.51   0.24 1.00     2036     1983
##    y1[2]    0.43   0.41 1.10 1.07  -1.37   2.26 1.00     1898     1803
##    y1[3]   -2.92  -2.93 1.15 1.11  -4.77  -1.01 1.00     1939     2002
##    y1[4]    5.71   5.74 1.12 1.08   3.88   7.56 1.00     1933     2014
##    y1[5]   -5.13  -5.12 1.16 1.13  -7.02  -3.21 1.00     1455     1332
##    y1[6]   -1.85  -1.86 1.14 1.17  -3.72  -0.04 1.00     1612     1781
##    y1[7]   -1.63  -1.66 1.10 1.06  -3.44   0.19 1.00     1929     1822
##    y1[8]    6.33   6.32 1.17 1.11   4.38   8.21 1.00     1615     1724
##    y1[9]   -0.20  -0.21 1.15 1.14  -2.07   1.67 1.00     1385     1956
##    y1[10]   3.96   3.97 1.10 1.04   2.13   5.81 1.00     1926     1919
##    y1[11]   3.28   3.26 1.13 1.08   1.43   5.18 1.00     2023     2097
##    y1[12]   2.74   2.72 1.10 1.06   0.96   4.54 1.00     1711     1853
##    y1[13]   7.94   7.92 1.16 1.12   6.15   9.86 1.00     2091     1889
##    y1[14]   2.70   2.72 1.23 1.16   0.60   4.68 1.00      914     1402
##    y1[15]  -2.72  -2.71 1.16 1.17  -4.62  -0.79 1.00     1476     1529
##    y1[16]   2.11   2.13 1.13 1.11   0.20   3.92 1.00     1711     1702
##    y1[17]   1.28   1.30 1.11 1.08  -0.60   3.07 1.00     1739     1915
##    y1[18]   6.78   6.81 1.12 1.10   4.96   8.64 1.00     1990     1961
##    y1[19]  -0.93  -0.93 1.07 1.03  -2.67   0.81 1.00     1955     1738
##    y1[20]   2.94   2.91 1.09 1.04   1.12   4.75 1.00     2102     1973
##    y1[21]   2.40   2.43 1.14 1.10   0.46   4.25 1.00     1619     1723
##    y1[22]   2.05   2.05 1.09 1.05   0.29   3.86 1.00     1978     1689
##    y1[23]  -5.63  -5.64 1.17 1.14  -7.50  -3.75 1.00     1832     1857
##    y1[24]  -0.12  -0.11 1.10 1.06  -1.91   1.71 1.00     1880     1930
##    y1[25]  -1.31  -1.33 1.12 1.07  -3.14   0.53 1.00     1530     1886
##    y1[26]   1.37   1.37 1.08 1.06  -0.37   3.16 1.00     2004     1931
##    y1[27]  -2.48  -2.48 1.11 1.07  -4.36  -0.66 1.00     1737     1932
##    y1[28]  -1.86  -1.85 1.18 1.15  -3.75   0.10 1.00      912     1227
##    y1[29]   6.78   6.75 1.14 1.08   4.89   8.68 1.00     1799     1733
##    y1[30]   4.89   4.89 1.14 1.10   3.05   6.76 1.00     1640     1839
##    m1[1]   -1.61  -1.61 0.27 0.27  -2.05  -1.14 1.00     1472     1389
##    m1[2]    0.45   0.45 0.31 0.31  -0.06   0.95 1.00      600      964
##    m1[3]   -2.92  -2.91 0.34 0.33  -3.49  -2.35 1.00     1152     1412
##    m1[4]    5.71   5.70 0.31 0.31   5.20   6.22 1.00     1683     1322
##    m1[5]   -5.14  -5.13 0.46 0.45  -5.91  -4.38 1.01      506      623
##    m1[6]   -1.86  -1.87 0.34 0.34  -2.44  -1.30 1.00      748     1297
##    m1[7]   -1.59  -1.59 0.25 0.24  -2.01  -1.18 1.00     1874     1488
##    m1[8]    6.31   6.32 0.41 0.39   5.63   6.98 1.00      615      648
##    m1[9]   -0.22  -0.23 0.37 0.38  -0.82   0.37 1.01      485      894
##    m1[10]   3.90   3.91 0.30 0.29   3.42   4.39 1.00      919     1246
##    m1[11]   3.28   3.28 0.23 0.22   2.92   3.65 1.00     1813     1508
##    m1[12]   2.74   2.74 0.25 0.24   2.33   3.15 1.01      570      454
##    m1[13]   7.93   7.92 0.41 0.40   7.27   8.61 1.00     1684     1267
##    m1[14]   2.67   2.69 0.54 0.54   1.80   3.53 1.01      271      279
##    m1[15]  -2.73  -2.72 0.35 0.34  -3.32  -2.15 1.01      432      422
##    m1[16]   2.06   2.04 0.41 0.41   1.39   2.70 1.01      403      701
##    m1[17]   1.28   1.28 0.28 0.27   0.83   1.72 1.01      358      344
##    m1[18]   6.82   6.80 0.36 0.35   6.25   7.40 1.00     1857     1269
##    m1[19]  -0.90  -0.91 0.23 0.22  -1.28  -0.53 1.00     1905     1365
##    m1[20]   2.91   2.91 0.24 0.23   2.52   3.30 1.00     1706     1260
##    m1[21]   2.36   2.37 0.34 0.34   1.80   2.90 1.00      499      985
##    m1[22]   2.05   2.06 0.21 0.20   1.71   2.41 1.00     1019      865
##    m1[23]  -5.62  -5.62 0.42 0.40  -6.31  -4.91 1.00     1216     1388
##    m1[24]  -0.12  -0.12 0.21 0.20  -0.47   0.23 1.01     1024      814
##    m1[25]  -1.31  -1.32 0.35 0.35  -1.89  -0.74 1.00      631     1087
##    m1[26]   1.37   1.37 0.20 0.20   1.04   1.71 1.00     2062     1444
##    m1[27]  -2.46  -2.46 0.31 0.30  -2.98  -1.94 1.01      556      603
##    m1[28]  -1.83  -1.81 0.51 0.51  -2.66  -1.01 1.01      277      277
##    m1[29]   6.79   6.79 0.37 0.37   6.17   7.40 1.00     1143     1161
##    m1[30]   4.91   4.92 0.38 0.37   4.29   5.54 1.01      443      450

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANOVA

tb=tibble(c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)

qplot(data=tb,c,y,geom='boxplot')

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -52.9384 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       14      -16.6315    0.00024107    0.00123899      0.9856      0.9856       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -16.63
##    b[1]      -0.18
##    b[2]       2.15
##    b[3]      -2.35
##    s          1.06
##    y1[1]     -0.42
##    y1[2]      3.05
##    y1[3]      0.80
##    y1[4]     -0.22
##    y1[5]     -1.48
##    y1[6]     -0.75
##    y1[7]      3.59
##    y1[8]      1.46
##    y1[9]      2.83
##    y1[10]    -0.17
##    y1[11]     0.86
##    y1[12]    -1.61
##    y1[13]    -0.23
##    y1[14]    -0.40
##    y1[15]     1.04
##    y1[16]     2.47
##    y1[17]     0.66
##    y1[18]     0.42
##    y1[19]    -3.97
##    y1[20]    -1.22
##    y1[21]    -1.89
##    y1[22]     0.18
##    y1[23]    -1.63
##    y1[24]    -3.74
##    y1[25]     1.34
##    y1[26]     2.07
##    y1[27]    -0.87
##    y1[28]    -3.20
##    y1[29]     0.99
##    y1[30]    -3.61
##    m1[1]     -2.53
##    m1[2]      1.97
##    m1[3]      1.97
##    m1[4]     -0.18
##    m1[5]     -0.18
##    m1[6]     -2.53
##    m1[7]      1.97
##    m1[8]      1.97
##    m1[9]      1.97
##    m1[10]    -2.53
##    m1[11]    -0.18
##    m1[12]    -2.53
##    m1[13]    -0.18
##    m1[14]    -2.53
##    m1[15]     1.97
##    m1[16]     1.97
##    m1[17]     1.97
##    m1[18]    -0.18
##    m1[19]    -2.53
##    m1[20]    -0.18
##    m1[21]    -0.18
##    m1[22]    -0.18
##    m1[23]    -2.53
##    m1[24]    -2.53
##    m1[25]     1.97
##    m1[26]     1.97
##    m1[27]    -2.53
##    m1[28]    -0.18
##    m1[29]    -0.18
##    m1[30]    -2.53
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -18.70 -18.36 1.52 1.34 -21.56 -16.91 1.00      629      943
##    b[1]    -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    b[2]     2.15   2.17 0.52 0.50   1.30   3.02 1.00      712      744
##    b[3]    -2.33  -2.34 0.52 0.52  -3.19  -1.48 1.01      753      807
##    s        1.17   1.15 0.17 0.16   0.93   1.48 1.00     1432     1442
##    y1[1]   -2.49  -2.45 1.24 1.25  -4.53  -0.50 1.00     1731     1821
##    y1[2]    1.98   1.97 1.26 1.21  -0.11   4.02 1.00     1894     2001
##    y1[3]    1.94   1.94 1.24 1.20  -0.07   3.95 1.00     1901     1923
##    y1[4]   -0.15  -0.15 1.25 1.19  -2.20   1.89 1.00     1731     1612
##    y1[5]   -0.21  -0.19 1.27 1.26  -2.33   1.80 1.00     1691     2101
##    y1[6]   -2.58  -2.57 1.24 1.21  -4.57  -0.54 1.00     1861     1842
##    y1[7]    1.97   1.98 1.27 1.26  -0.09   4.11 1.00     1800     1806
##    y1[8]    1.94   1.90 1.22 1.22  -0.02   3.99 1.00     1795     1854
##    y1[9]    1.92   1.94 1.21 1.15  -0.03   3.97 1.00     2008     1875
##    y1[10]  -2.53  -2.53 1.24 1.21  -4.54  -0.53 1.00     1832     1637
##    y1[11]  -0.20  -0.19 1.22 1.19  -2.21   1.84 1.00     1823     1841
##    y1[12]  -2.51  -2.52 1.25 1.26  -4.59  -0.48 1.00     2305     1876
##    y1[13]  -0.20  -0.19 1.23 1.25  -2.22   1.81 1.00     1622     1875
##    y1[14]  -2.56  -2.55 1.26 1.22  -4.64  -0.49 1.00     2057     1615
##    y1[15]   1.99   1.99 1.23 1.20   0.00   4.03 1.00     2121     1841
##    y1[16]   1.96   1.96 1.22 1.19  -0.06   3.92 1.00     1903     1990
##    y1[17]   1.98   1.99 1.24 1.22  -0.05   4.02 1.00     1828     1984
##    y1[18]  -0.18  -0.20 1.22 1.23  -2.12   1.81 1.00     1803     1916
##    y1[19]  -2.53  -2.49 1.26 1.20  -4.59  -0.48 1.00     1830     1915
##    y1[20]  -0.18  -0.20 1.27 1.23  -2.18   1.90 1.00     1640     1390
##    y1[21]  -0.18  -0.19 1.26 1.23  -2.19   1.89 1.00     1918     1348
##    y1[22]  -0.17  -0.20 1.28 1.21  -2.37   1.89 1.00     1754     1819
##    y1[23]  -2.55  -2.54 1.24 1.21  -4.57  -0.62 1.00     1982     1713
##    y1[24]  -2.51  -2.51 1.24 1.19  -4.56  -0.51 1.00     1898     1796
##    y1[25]   1.95   1.98 1.25 1.22  -0.01   4.02 1.00     1787     1886
##    y1[26]   1.99   1.98 1.26 1.18  -0.06   4.04 1.00     2025     1698
##    y1[27]  -2.52  -2.55 1.25 1.20  -4.50  -0.45 1.00     1854     1934
##    y1[28]  -0.22  -0.23 1.24 1.18  -2.27   1.82 1.00     1973     1880
##    y1[29]  -0.20  -0.19 1.26 1.23  -2.28   1.83 1.00     1932     1931
##    y1[30]  -2.50  -2.48 1.22 1.18  -4.54  -0.54 1.00     1937     1821
##    m1[1]   -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[2]    1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[3]    1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[4]   -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[5]   -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[6]   -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[7]    1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[8]    1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[9]    1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[10]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[11]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[12]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[13]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[14]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[15]   1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[16]   1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[17]   1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[18]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[19]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[20]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[21]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[22]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[23]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[24]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[25]   1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[26]   1.96   1.97 0.37 0.36   1.36   2.57 1.00     1130     1219
##    m1[27]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346
##    m1[28]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[29]  -0.19  -0.20 0.37 0.35  -0.78   0.41 1.00      820     1042
##    m1[30]  -2.52  -2.53 0.37 0.36  -3.14  -1.89 1.01     1121     1346

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANCOVA

tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))

f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)

lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])

qplot(data=tb,x,y,shape=c,size=I(2))

plot(tb$x,tb$y,col=1+lv[factor(tb$c)])

qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -451.38 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       43      -9.28787   6.39702e-05   0.000959384           1           1       49    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -9.29
##    b[1]       0.40
##    b[2]       0.93
##    b[3]       1.62
##    b[4]      -2.19
##    s          0.83
##    y1[1]      1.90
##    y1[2]      5.40
##    y1[3]      2.89
##    y1[4]     -1.32
##    y1[5]      0.37
##    y1[6]      7.28
##    y1[7]     -1.70
##    y1[8]      9.26
##    y1[9]      0.85
##    y1[10]     7.72
##    y1[11]     4.09
##    y1[12]     8.41
##    y1[13]     8.72
##    y1[14]     3.41
##    y1[15]     7.33
##    y1[16]     6.15
##    y1[17]     6.99
##    y1[18]     1.88
##    y1[19]     3.16
##    y1[20]     3.86
##    y1[21]     4.37
##    y1[22]     2.42
##    y1[23]     7.04
##    y1[24]     3.61
##    y1[25]     4.29
##    y1[26]    -0.75
##    y1[27]     0.90
##    y1[28]     6.71
##    y1[29]     6.05
##    y1[30]     8.05
##    m1[1]      3.51
##    m1[2]      4.76
##    m1[3]      3.82
##    m1[4]     -0.35
##    m1[5]      0.86
##    m1[6]      7.81
##    m1[7]     -1.12
##    m1[8]      7.89
##    m1[9]     -0.10
##    m1[10]     6.65
##    m1[11]     3.19
##    m1[12]     8.32
##    m1[13]     9.06
##    m1[14]     5.43
##    m1[15]     6.77
##    m1[16]     5.75
##    m1[17]     7.63
##    m1[18]     2.18
##    m1[19]     2.84
##    m1[20]     4.01
##    m1[21]     5.32
##    m1[22]     2.75
##    m1[23]     6.85
##    m1[24]     4.09
##    m1[25]     5.75
##    m1[26]     0.37
##    m1[27]     1.99
##    m1[28]     5.94
##    m1[29]     6.90
##    m1[30]     6.51
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.08 -11.75 1.65 1.44 -15.28 -10.05 1.00      622      840
##    b[1]     0.38   0.37 0.44 0.43  -0.31   1.12 1.01      607      702
##    b[2]     0.93   0.93 0.07 0.07   0.82   1.04 1.00     1148     1151
##    b[3]     1.63   1.62 0.40 0.38   0.98   2.29 1.01      683      961
##    b[4]    -2.20  -2.19 0.40 0.39  -2.88  -1.55 1.00      643      767
##    s        0.93   0.92 0.14 0.13   0.74   1.16 1.00     1868     1480
##    y1[1]    3.49   3.47 1.00 1.01   1.87   5.14 1.00     1762     2074
##    y1[2]    4.76   4.74 0.97 0.92   3.24   6.38 1.00     1845     2014
##    y1[3]    3.80   3.81 1.01 0.99   2.19   5.46 1.00     2019     1922
##    y1[4]   -0.35  -0.37 1.01 0.98  -2.00   1.30 1.00     1975     1971
##    y1[5]    0.89   0.90 1.04 1.04  -0.85   2.56 1.01     1270     1666
##    y1[6]    7.80   7.79 0.97 0.98   6.23   9.41 1.00     1878     1846
##    y1[7]   -1.14  -1.14 1.04 1.04  -2.88   0.52 1.00     1945     1788
##    y1[8]    7.89   7.87 1.01 0.98   6.22   9.57 1.00     1792     1917
##    y1[9]   -0.16  -0.16 0.99 0.95  -1.81   1.46 1.00     1830     1911
##    y1[10]   6.67   6.68 1.00 0.96   5.03   8.31 1.00     2054     1959
##    y1[11]   3.18   3.19 0.99 0.97   1.55   4.84 1.00     2005     1857
##    y1[12]   8.33   8.32 1.01 1.01   6.72  10.00 1.00     2094     1841
##    y1[13]   9.09   9.12 1.00 0.95   7.38  10.72 1.00     1775     1559
##    y1[14]   5.44   5.44 0.97 1.00   3.94   7.03 1.00     1713     2052
##    y1[15]   6.79   6.80 1.01 0.97   5.16   8.42 1.00     2073     1972
##    y1[16]   5.76   5.81 0.98 0.96   4.15   7.33 1.00     1994     1930
##    y1[17]   7.66   7.64 1.00 0.95   6.07   9.33 1.00     1869     1736
##    y1[18]   2.18   2.16 1.04 1.01   0.54   3.91 1.00     1843     1630
##    y1[19]   2.85   2.83 0.99 0.93   1.27   4.47 1.00     2051     1802
##    y1[20]   4.01   3.98 1.01 0.98   2.44   5.67 1.00     1743     1618
##    y1[21]   5.31   5.31 0.99 0.95   3.67   6.91 1.00     1851     1858
##    y1[22]   2.71   2.71 0.98 0.99   1.12   4.32 1.00     1679     1879
##    y1[23]   6.86   6.83 0.99 0.96   5.18   8.51 1.00     2045     1800
##    y1[24]   4.05   4.03 1.00 0.98   2.46   5.68 1.00     1981     2014
##    y1[25]   5.75   5.77 0.99 1.00   4.12   7.33 1.00     1930     1721
##    y1[26]   0.37   0.35 1.00 0.98  -1.23   2.04 1.00     1969     1866
##    y1[27]   1.99   1.98 0.98 0.99   0.43   3.58 1.00     1910     1929
##    y1[28]   5.93   5.95 0.96 0.92   4.31   7.45 1.00     1921     1751
##    y1[29]   6.90   6.89 0.98 0.97   5.29   8.48 1.00     1988     1892
##    y1[30]   6.50   6.50 1.02 1.04   4.78   8.15 1.00     1544     1706
##    m1[1]    3.50   3.50 0.30 0.29   3.00   3.99 1.00     1251     1444
##    m1[2]    4.76   4.75 0.29 0.28   4.31   5.25 1.00      684      761
##    m1[3]    3.81   3.82 0.31 0.30   3.31   4.31 1.00     1232     1489
##    m1[4]   -0.37  -0.38 0.36 0.35  -0.95   0.20 1.00     1324     1068
##    m1[5]    0.84   0.83 0.42 0.40   0.19   1.55 1.01      595      662
##    m1[6]    7.81   7.80 0.35 0.34   7.25   8.39 1.00     1415     1426
##    m1[7]   -1.14  -1.15 0.39 0.38  -1.77  -0.51 1.00     1308     1014
##    m1[8]    7.90   7.91 0.33 0.32   7.36   8.42 1.01     1266     1436
##    m1[9]   -0.12  -0.13 0.35 0.33  -0.68   0.44 1.00     1333     1149
##    m1[10]   6.66   6.66 0.30 0.29   6.17   7.15 1.00     1367     1485
##    m1[11]   3.19   3.19 0.35 0.33   2.64   3.78 1.00     1383     1294
##    m1[12]   8.32   8.32 0.37 0.36   7.73   8.93 1.00     1550     1458
##    m1[13]   9.07   9.08 0.37 0.36   8.45   9.67 1.01     1200     1429
##    m1[14]   5.42   5.42 0.29 0.28   4.97   5.91 1.00      781      894
##    m1[15]   6.77   6.78 0.30 0.29   6.28   7.27 1.00     1358     1386
##    m1[16]   5.75   5.74 0.29 0.28   5.27   6.23 1.00     1402     1453
##    m1[17]   7.64   7.65 0.32 0.31   7.11   8.16 1.01     1285     1448
##    m1[18]   2.17   2.17 0.39 0.38   1.54   2.84 1.00     1354     1196
##    m1[19]   2.83   2.83 0.36 0.35   2.26   3.44 1.00     1376     1181
##    m1[20]   4.00   3.99 0.30 0.29   3.54   4.51 1.01      615      706
##    m1[21]   5.31   5.32 0.36 0.35   4.71   5.90 1.00     1193     1401
##    m1[22]   2.74   2.73 0.33 0.32   2.21   3.30 1.01      575      674
##    m1[23]   6.85   6.84 0.31 0.31   6.35   7.37 1.00     1143     1286
##    m1[24]   4.08   4.09 0.32 0.31   3.57   4.59 1.00     1223     1469
##    m1[25]   5.75   5.75 0.29 0.28   5.28   6.24 1.00     1402     1453
##    m1[26]   0.34   0.33 0.33 0.32  -0.20   0.89 1.00     1344     1164
##    m1[27]   1.97   1.98 0.30 0.29   1.49   2.45 1.00     1318     1243
##    m1[28]   5.93   5.93 0.29 0.29   5.47   6.43 1.00      884      946
##    m1[29]   6.90   6.89 0.31 0.31   6.40   7.42 1.00     1157     1360
##    m1[30]   6.50   6.52 0.41 0.40   5.81   7.18 1.00     1163     1264

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
  geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

interaction of variable

n=50
mdl=cmdstan_model('./ex5-1.stan')

tb=tibble(x=runif(n,-3,3),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,x,y,col=ca),
             qplot(data=tb,x,y,col=cb),ncol=2)

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)

fn(f0)
## Initial log joint probability = -119.88 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -27.5817   0.000129179   0.000855749           1           1       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -29.64 -29.28 1.52 1.28 -32.48 -27.88 1.01      908     1002
##    b[1]     0.71   0.71 0.22 0.21   0.34   1.09 1.00      586      987
##    b[2]     0.97   0.97 0.09 0.09   0.81   1.12 1.00     2127     1383
##    b[3]     0.51   0.52 0.32 0.33  -0.01   1.03 1.00      515      861
##    s        1.11   1.11 0.12 0.12   0.93   1.33 1.00     1529     1448
##    y1[1]   -1.71  -1.75 1.17 1.15  -3.59   0.27 1.00     1932     2031
##    y1[2]   -1.79  -1.78 1.16 1.08  -3.71   0.13 1.00     1943     1966
##    y1[3]    0.49   0.51 1.16 1.18  -1.44   2.34 1.00     1921     1933
##    y1[4]    1.18   1.16 1.12 1.09  -0.64   3.06 1.00     2034     1929
##    y1[5]    1.19   1.20 1.14 1.10  -0.72   3.02 1.00     2003     1896
##    y1[6]   -0.85  -0.86 1.16 1.11  -2.82   1.09 1.00     1912     1954
##    y1[7]    1.02   1.01 1.12 1.11  -0.79   2.87 1.00     1765     1909
##    y1[8]    1.24   1.24 1.14 1.09  -0.62   3.14 1.00     2074     1958
##    y1[9]    3.20   3.15 1.15 1.17   1.34   5.14 1.00     1808     1850
##    y1[10]   2.52   2.54 1.14 1.12   0.60   4.40 1.00     1901     1928
##    y1[11]   3.44   3.47 1.17 1.12   1.42   5.32 1.00     2154     1972
##    y1[12]   2.97   3.00 1.19 1.21   0.99   4.87 1.00     1419     1862
##    y1[13]  -1.15  -1.13 1.15 1.15  -3.12   0.67 1.00     2113     2096
##    y1[14]   2.75   2.73 1.17 1.12   0.84   4.66 1.00     1741     1785
##    y1[15]   3.46   3.47 1.19 1.18   1.50   5.39 1.00     2056     1894
##    y1[16]  -1.77  -1.78 1.13 1.10  -3.60   0.15 1.00     1818     1734
##    y1[17]  -1.00  -0.99 1.13 1.07  -2.86   0.79 1.00     1922     1823
##    y1[18]   1.58   1.57 1.16 1.17  -0.30   3.45 1.00     2202     2020
##    y1[19]   0.37   0.36 1.17 1.13  -1.57   2.34 1.00     1797     1908
##    y1[20]   0.59   0.59 1.17 1.19  -1.36   2.45 1.00     1830     1763
##    y1[21]   2.92   2.88 1.15 1.16   1.09   4.85 1.00     2047     1966
##    y1[22]  -1.38  -1.37 1.15 1.15  -3.21   0.51 1.00     1734     1811
##    y1[23]   0.78   0.79 1.13 1.10  -1.10   2.60 1.00     1638     1871
##    y1[24]  -2.03  -2.02 1.16 1.18  -3.93  -0.18 1.00     1958     1859
##    y1[25]  -1.47  -1.44 1.14 1.16  -3.38   0.41 1.00     2068     1934
##    y1[26]   3.38   3.40 1.18 1.19   1.38   5.33 1.00     1901     1835
##    y1[27]   0.40   0.40 1.14 1.10  -1.47   2.27 1.00     1832     1932
##    y1[28]   1.78   1.78 1.12 1.11  -0.06   3.60 1.00     1997     2024
##    y1[29]  -0.67  -0.67 1.16 1.14  -2.55   1.23 1.00     1801     1882
##    y1[30]  -1.75  -1.75 1.16 1.14  -3.69   0.14 1.00     1924     2004
##    y1[31]   0.99   1.00 1.13 1.10  -0.91   2.90 1.00     2056     1991
##    y1[32]   0.10   0.12 1.15 1.11  -1.83   1.95 1.00     2013     1839
##    y1[33]   0.59   0.61 1.15 1.16  -1.26   2.46 1.00     1846     1807
##    y1[34]  -1.33  -1.37 1.19 1.18  -3.26   0.68 1.00     1918     1929
##    y1[35]   3.23   3.21 1.16 1.13   1.37   5.18 1.00     1996     1886
##    y1[36]   2.97   2.96 1.15 1.12   1.16   4.93 1.00     1925     1762
##    y1[37]  -1.10  -1.12 1.20 1.21  -3.05   0.95 1.00     1940     1968
##    y1[38]   0.14   0.13 1.16 1.14  -1.71   2.06 1.00     1879     1947
##    y1[39]   0.56   0.59 1.16 1.16  -1.43   2.48 1.00     2077     1884
##    y1[40]   0.16   0.14 1.12 1.11  -1.63   2.00 1.00     1842     1972
##    y1[41]   1.19   1.20 1.14 1.10  -0.72   3.01 1.00     1972     1878
##    y1[42]   2.63   2.63 1.14 1.14   0.75   4.54 1.00     1985     2012
##    y1[43]  -1.29  -1.29 1.13 1.13  -3.10   0.55 1.00     2029     2053
##    y1[44]   0.92   0.93 1.14 1.12  -0.98   2.70 1.00     1961     1970
##    y1[45]  -1.20  -1.20 1.15 1.12  -3.07   0.69 1.00     1749     1754
##    y1[46]   1.16   1.19 1.15 1.13  -0.74   2.97 1.00     1931     1921
##    y1[47]  -2.15  -2.18 1.17 1.15  -4.03  -0.26 1.00     1829     2034
##    y1[48]   3.41   3.41 1.17 1.13   1.50   5.35 1.00     1860     1959
##    y1[49]   3.49   3.51 1.22 1.19   1.38   5.44 1.00     1765     1886
##    y1[50]   1.17   1.19 1.15 1.14  -0.74   3.05 1.00     1920     2010
##    m1[1]   -1.69  -1.69 0.28 0.27  -2.15  -1.23 1.00     1455     1756
##    m1[2]   -1.83  -1.83 0.29 0.28  -2.30  -1.35 1.00     1523     1757
##    m1[3]    0.48   0.48 0.24 0.23   0.08   0.87 1.00      987     1221
##    m1[4]    1.15   1.15 0.23 0.22   0.76   1.52 1.00     1050     1091
##    m1[5]    1.16   1.16 0.23 0.23   0.77   1.55 1.00      564      947
##    m1[6]   -0.83  -0.83 0.24 0.24  -1.23  -0.44 1.00      991     1298
##    m1[7]    1.05   1.05 0.23 0.22   0.67   1.42 1.00     1035     1104
##    m1[8]    1.24   1.24 0.23 0.22   0.86   1.62 1.00     1060     1110
##    m1[9]    3.19   3.19 0.30 0.28   2.69   3.68 1.00     1704     1480
##    m1[10]   2.58   2.58 0.26 0.25   2.13   3.01 1.00     1466     1308
##    m1[11]   3.46   3.47 0.31 0.30   2.93   3.98 1.00     1808     1503
##    m1[12]   3.00   3.01 0.35 0.34   2.43   3.56 1.00      706     1023
##    m1[13]  -1.13  -1.14 0.25 0.25  -1.55  -0.72 1.00     1172     1592
##    m1[14]   2.73   2.74 0.33 0.32   2.19   3.27 1.00      677     1043
##    m1[15]   3.40   3.40 0.31 0.30   2.87   3.91 1.00     1784     1529
##    m1[16]  -1.81  -1.81 0.29 0.28  -2.28  -1.33 1.00     1517     1757
##    m1[17]  -1.01  -1.01 0.25 0.24  -1.41  -0.61 1.00     1091     1539
##    m1[18]   1.59   1.59 0.23 0.22   1.21   1.97 1.00     1129     1134
##    m1[19]   0.33   0.33 0.24 0.24  -0.07   0.72 1.00      983     1322
##    m1[20]   0.60   0.60 0.22 0.21   0.23   0.97 1.00      601     1008
##    m1[21]   2.96   2.96 0.28 0.27   2.48   3.42 1.00     1610     1414
##    m1[22]  -1.36  -1.36 0.34 0.33  -1.91  -0.82 1.00     1116     1530
##    m1[23]   0.74   0.74 0.22 0.21   0.37   1.12 1.00      583      987
##    m1[24]  -2.00  -2.00 0.30 0.29  -2.49  -1.50 1.00     1607     1663
##    m1[25]  -1.47  -1.47 0.27 0.26  -1.91  -1.03 1.00     1346     1621
##    m1[26]   3.35   3.36 0.37 0.37   2.73   3.96 1.00      739      911
##    m1[27]   0.35   0.35 0.24 0.24  -0.06   0.73 1.00      984     1322
##    m1[28]   1.79   1.79 0.23 0.23   1.41   2.18 1.00     1186     1175
##    m1[29]  -0.65  -0.65 0.29 0.28  -1.14  -0.20 1.00     1024     1555
##    m1[30]  -1.75  -1.75 0.29 0.27  -2.22  -1.28 1.00     1487     1676
##    m1[31]   1.03   1.04 0.23 0.22   0.65   1.41 1.00     1033     1103
##    m1[32]   0.14   0.14 0.25 0.24  -0.28   0.54 1.00      981     1350
##    m1[33]   0.62   0.62 0.22 0.22   0.25   0.99 1.00      598     1008
##    m1[34]  -1.31  -1.31 0.33 0.32  -1.85  -0.77 1.00     1108     1482
##    m1[35]   3.22   3.22 0.30 0.28   2.71   3.71 1.00     1717     1414
##    m1[36]   2.96   2.96 0.28 0.27   2.48   3.42 1.00     1609     1414
##    m1[37]  -1.14  -1.14 0.32 0.31  -1.66  -0.62 1.00     1086     1512
##    m1[38]   0.13   0.13 0.22 0.22  -0.23   0.49 1.00      681     1017
##    m1[39]   0.55   0.55 0.24 0.23   0.16   0.93 1.00      992     1221
##    m1[40]   0.13   0.13 0.22 0.22  -0.23   0.49 1.00      681     1017
##    m1[41]   1.23   1.23 0.24 0.24   0.84   1.62 1.00      564      834
##    m1[42]   2.64   2.64 0.27 0.26   2.19   3.08 1.00     1490     1399
##    m1[43]  -1.29  -1.29 0.26 0.25  -1.72  -0.86 1.00     1259     1570
##    m1[44]   0.91   0.91 0.23 0.22   0.53   1.29 1.00      572      989
##    m1[45]  -1.20  -1.20 0.25 0.25  -1.62  -0.78 1.00     1217     1445
##    m1[46]   1.19   1.19 0.23 0.22   0.81   1.57 1.00     1054     1117
##    m1[47]  -2.13  -2.12 0.31 0.30  -2.64  -1.62 1.00     1671     1634
##    m1[48]   3.43   3.43 0.38 0.37   2.80   4.04 1.00      748      907
##    m1[49]   3.45   3.45 0.38 0.37   2.82   4.06 1.00      750      907
##    m1[50]   1.18   1.18 0.23 0.22   0.80   1.56 1.00     1053     1117

fn(f1)
## Initial log joint probability = -1090.15 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -26.7666   0.000133064     0.0011364           1           1       20    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -29.37 -29.02 1.67 1.42 -32.63 -27.37 1.00      866     1174
##    b[1]     0.47   0.46 0.30 0.30  -0.03   0.95 1.01      529      738
##    b[2]     0.98   0.98 0.09 0.09   0.83   1.13 1.00     2201     1535
##    b[3]     0.50   0.52 0.33 0.32  -0.06   1.04 1.00      675      844
##    b[4]     0.40   0.41 0.31 0.32  -0.11   0.90 1.01      695      816
##    s        1.11   1.10 0.12 0.12   0.93   1.33 1.00     2074     1631
##    y1[1]   -1.54  -1.53 1.17 1.13  -3.44   0.35 1.00     1817     1911
##    y1[2]   -1.67  -1.71 1.13 1.09  -3.43   0.21 1.00     2168     1897
##    y1[3]    0.23   0.26 1.16 1.14  -1.65   2.13 1.00     1931     1933
##    y1[4]    1.30   1.29 1.14 1.15  -0.55   3.18 1.00     2087     1927
##    y1[5]    1.33   1.34 1.16 1.14  -0.58   3.20 1.00     1757     1669
##    y1[6]   -0.69  -0.70 1.17 1.17  -2.62   1.25 1.00     1925     1910
##    y1[7]    1.21   1.19 1.12 1.11  -0.60   2.99 1.00     1734     1923
##    y1[8]    0.97   0.99 1.17 1.15  -0.99   2.93 1.00     1611     1928
##    y1[9]    3.31   3.28 1.17 1.12   1.38   5.21 1.00     1884     1791
##    y1[10]   2.74   2.74 1.17 1.11   0.86   4.69 1.00     1847     1825
##    y1[11]   3.62   3.61 1.16 1.13   1.67   5.57 1.00     1980     1917
##    y1[12]   3.24   3.27 1.16 1.15   1.32   5.16 1.00     1857     1833
##    y1[13]  -0.97  -1.01 1.17 1.13  -2.89   0.95 1.00     2076     2150
##    y1[14]   2.49   2.51 1.17 1.13   0.54   4.38 1.00     1913     2011
##    y1[15]   3.19   3.22 1.21 1.18   1.15   5.12 1.00     1884     1971
##    y1[16]  -1.66  -1.67 1.18 1.22  -3.61   0.28 1.00     1997     1927
##    y1[17]  -0.87  -0.85 1.16 1.16  -2.80   1.03 1.00     2010     2012
##    y1[18]   1.37   1.38 1.16 1.19  -0.52   3.31 1.00     2067     1973
##    y1[19]   0.01   0.05 1.15 1.14  -1.88   1.85 1.00     1700     1769
##    y1[20]   0.79   0.77 1.13 1.11  -1.06   2.68 1.00     1715     1933
##    y1[21]   2.69   2.69 1.16 1.18   0.79   4.61 1.00     2086     1821
##    y1[22]  -1.62  -1.62 1.20 1.18  -3.60   0.40 1.00     1995     1803
##    y1[23]   0.48   0.49 1.16 1.13  -1.47   2.38 1.01     1547     1858
##    y1[24]  -2.22  -2.23 1.19 1.15  -4.17  -0.29 1.00     1727     1747
##    y1[25]  -1.32  -1.31 1.19 1.17  -3.34   0.57 1.00     1910     2016
##    y1[26]   3.11   3.11 1.19 1.18   1.19   5.09 1.00     1728     1617
##    y1[27]   0.51   0.50 1.16 1.11  -1.41   2.36 1.00     1914     1932
##    y1[28]   1.98   1.97 1.15 1.14   0.09   3.88 1.00     1921     1855
##    y1[29]  -0.49  -0.50 1.19 1.19  -2.43   1.52 1.00     1859     1639
##    y1[30]  -2.05  -2.03 1.17 1.15  -3.95  -0.11 1.00     1969     2045
##    y1[31]   1.16   1.15 1.16 1.14  -0.72   3.04 1.00     1499     1636
##    y1[32]   0.28   0.27 1.14 1.10  -1.63   2.13 1.00     1706     1788
##    y1[33]   0.74   0.73 1.13 1.12  -1.04   2.63 1.00     1732     1805
##    y1[34]  -1.22  -1.19 1.14 1.09  -3.10   0.64 1.00     2087     1924
##    y1[35]   3.37   3.38 1.15 1.12   1.51   5.23 1.00     1947     1853
##    y1[36]   3.09   3.09 1.12 1.12   1.23   4.86 1.00     2113     2102
##    y1[37]  -1.04  -1.03 1.17 1.14  -3.04   0.85 1.00     1773     1886
##    y1[38]   0.25   0.24 1.15 1.11  -1.59   2.09 1.00     1941     1890
##    y1[39]   0.65   0.67 1.13 1.09  -1.19   2.51 1.00     1929     1820
##    y1[40]   0.28   0.26 1.14 1.17  -1.61   2.13 1.00     2086     2052
##    y1[41]   1.37   1.36 1.12 1.07  -0.45   3.31 1.00     1649     1801
##    y1[42]   2.42   2.42 1.14 1.17   0.54   4.27 1.00     1997     1654
##    y1[43]  -1.54  -1.51 1.16 1.14  -3.51   0.30 1.00     1586     1827
##    y1[44]   0.66   0.65 1.13 1.14  -1.18   2.52 1.00     1665     1723
##    y1[45]  -1.08  -1.06 1.15 1.12  -3.02   0.81 1.00     1975     2035
##    y1[46]   1.34   1.34 1.15 1.13  -0.57   3.25 1.00     1901     2058
##    y1[47]  -2.41  -2.39 1.19 1.11  -4.41  -0.42 1.00     1815     1752
##    y1[48]   3.19   3.16 1.16 1.12   1.32   5.09 1.00     1837     1836
##    y1[49]   3.66   3.64 1.20 1.20   1.74   5.66 1.00     1911     1892
##    y1[50]   0.90   0.88 1.15 1.15  -0.97   2.83 1.00     1804     2016
##    m1[1]   -1.56  -1.56 0.31 0.30  -2.05  -1.05 1.00     1968     1455
##    m1[2]   -1.70  -1.70 0.32 0.30  -2.21  -1.18 1.00     2025     1428
##    m1[3]    0.22   0.21 0.31 0.30  -0.30   0.72 1.00      959     1209
##    m1[4]    1.29   1.30 0.26 0.26   0.86   1.71 1.00     1266     1305
##    m1[5]    1.32   1.33 0.28 0.28   0.86   1.77 1.00     1022     1184
##    m1[6]   -0.69  -0.69 0.27 0.26  -1.13  -0.24 1.00     1613     1522
##    m1[7]    1.19   1.20 0.26 0.26   0.77   1.61 1.00     1253     1361
##    m1[8]    0.98   0.98 0.30 0.29   0.49   1.47 1.00      956     1156
##    m1[9]    3.36   3.36 0.32 0.30   2.83   3.90 1.00     1757     1734
##    m1[10]   2.74   2.74 0.29 0.28   2.25   3.22 1.00     1573     1484
##    m1[11]   3.63   3.63 0.33 0.33   3.09   4.21 1.00     1829     1587
##    m1[12]   3.19   3.18 0.38 0.38   2.57   3.81 1.00     1217     1554
##    m1[13]  -1.00  -1.00 0.28 0.28  -1.45  -0.53 1.00     1739     1524
##    m1[14]   2.51   2.51 0.38 0.38   1.90   3.14 1.00      687      881
##    m1[15]   3.16   3.16 0.35 0.34   2.58   3.74 1.00     1291     1468
##    m1[16]  -1.68  -1.68 0.32 0.30  -2.19  -1.16 1.00     2018     1407
##    m1[17]  -0.87  -0.87 0.28 0.27  -1.32  -0.41 1.00     1684     1610
##    m1[18]   1.34   1.34 0.30 0.29   0.86   1.83 1.00      973     1124
##    m1[19]   0.07   0.06 0.31 0.31  -0.45   0.58 1.00      964     1176
##    m1[20]   0.75   0.76 0.26 0.26   0.31   1.17 1.00     1085     1239
##    m1[21]   2.72   2.72 0.33 0.32   2.18   3.26 1.00     1181     1488
##    m1[22]  -1.65  -1.65 0.40 0.39  -2.31  -1.00 1.01     1144     1541
##    m1[23]   0.50   0.49 0.30 0.30   0.00   0.99 1.01      529      738
##    m1[24]  -2.27  -2.28 0.37 0.38  -2.86  -1.69 1.01      863     1073
##    m1[25]  -1.34  -1.34 0.30 0.29  -1.81  -0.83 1.00     1878     1436
##    m1[26]   3.14   3.13 0.42 0.41   2.45   3.84 1.00      773      889
##    m1[27]   0.48   0.49 0.27 0.27   0.03   0.92 1.00     1225     1256
##    m1[28]   1.95   1.94 0.26 0.26   1.50   2.38 1.00     1368     1366
##    m1[29]  -0.53  -0.52 0.32 0.31  -1.05  -0.02 1.00     1286     1306
##    m1[30]  -2.03  -2.03 0.36 0.36  -2.59  -1.45 1.01      808     1088
##    m1[31]   1.18   1.18 0.26 0.26   0.75   1.60 1.00     1252     1307
##    m1[32]   0.27   0.28 0.28 0.28  -0.19   0.72 1.00     1230     1270
##    m1[33]   0.77   0.78 0.26 0.27   0.33   1.19 1.00     1081     1256
##    m1[34]  -1.19  -1.19 0.36 0.35  -1.79  -0.61 1.00     1362     1290
##    m1[35]   3.39   3.39 0.32 0.31   2.85   3.93 1.00     1766     1734
##    m1[36]   3.13   3.12 0.30 0.29   2.61   3.64 1.00     1687     1700
##    m1[37]  -1.02  -1.01 0.35 0.34  -1.59  -0.46 1.00     1338     1272
##    m1[38]   0.28   0.29 0.26 0.26  -0.14   0.70 1.00     1197     1365
##    m1[39]   0.69   0.70 0.27 0.26   0.24   1.12 1.00     1222     1282
##    m1[40]   0.28   0.29 0.26 0.26  -0.14   0.70 1.00     1197     1365
##    m1[41]   1.39   1.40 0.28 0.28   0.92   1.84 1.00     1021     1194
##    m1[42]   2.40   2.40 0.32 0.30   1.89   2.92 1.00     1115     1541
##    m1[43]  -1.56  -1.56 0.33 0.33  -2.10  -1.01 1.01      715      951
##    m1[44]   0.67   0.66 0.31 0.30   0.17   1.16 1.01      532      709
##    m1[45]  -1.06  -1.06 0.29 0.28  -1.52  -0.58 1.00     1766     1524
##    m1[46]   1.33   1.34 0.26 0.25   0.90   1.75 1.00     1272     1305
##    m1[47]  -2.40  -2.40 0.38 0.38  -3.01  -1.79 1.01      894     1088
##    m1[48]   3.22   3.21 0.43 0.42   2.52   3.92 1.00      785      889
##    m1[49]   3.64   3.64 0.41 0.41   2.98   4.32 1.00     1283     1610
##    m1[50]   0.93   0.93 0.30 0.29   0.44   1.41 1.00      955     1104

fn(f2)
## Initial log joint probability = -81.3509 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       29      -25.1811   0.000169817     0.0004775      0.9558      0.9558       36    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -28.35 -27.98 1.87 1.60 -32.07 -26.03 1.00      697      849
##    b[1]     0.16   0.15 0.37 0.36  -0.43   0.74 1.00      539      828
##    b[2]     0.99   0.99 0.09 0.09   0.84   1.14 1.00     2080     1643
##    b[3]     1.16   1.17 0.52 0.52   0.26   1.99 1.00      454      691
##    b[4]     0.90   0.91 0.45 0.45   0.15   1.65 1.00      479      899
##    b[5]    -1.03  -1.04 0.65 0.66  -2.10   0.06 1.01      412      659
##    s        1.09   1.08 0.12 0.11   0.92   1.30 1.00     1812     1442
##    y1[1]   -1.37  -1.36 1.15 1.16  -3.24   0.50 1.00     2018     1968
##    y1[2]   -1.52  -1.51 1.16 1.14  -3.44   0.35 1.00     1948     1957
##    y1[3]    0.54   0.54 1.14 1.14  -1.39   2.42 1.00     1859     1932
##    y1[4]    1.08   1.10 1.15 1.11  -0.87   2.96 1.00     2028     2016
##    y1[5]    1.47   1.44 1.14 1.13  -0.40   3.32 1.00     1930     1870
##    y1[6]   -0.56  -0.57 1.09 1.08  -2.35   1.21 1.00     2019     1971
##    y1[7]    0.98   1.00 1.15 1.12  -0.94   2.79 1.00     2025     1915
##    y1[8]    1.32   1.33 1.16 1.13  -0.58   3.29 1.00     1735     1822
##    y1[9]    3.19   3.22 1.14 1.14   1.31   5.02 1.00     2010     1996
##    y1[10]   2.65   2.66 1.12 1.08   0.85   4.48 1.00     2023     1932
##    y1[11]   3.50   3.52 1.14 1.13   1.68   5.35 1.00     1860     1596
##    y1[12]   3.39   3.39 1.14 1.11   1.52   5.27 1.00     1764     1647
##    y1[13]  -0.86  -0.85 1.12 1.08  -2.72   1.03 1.00     2023     1857
##    y1[14]   2.19   2.20 1.15 1.12   0.32   4.12 1.00     1308     1768
##    y1[15]   3.53   3.55 1.17 1.17   1.64   5.44 1.00     1779     1850
##    y1[16]  -1.53  -1.52 1.16 1.14  -3.39   0.39 1.00     1927     1857
##    y1[17]  -0.73  -0.70 1.12 1.11  -2.63   1.10 1.00     1965     1848
##    y1[18]   1.68   1.64 1.17 1.11  -0.29   3.61 1.00     1764     1837
##    y1[19]   0.38   0.40 1.16 1.12  -1.59   2.21 1.00     1798     1823
##    y1[20]   0.89   0.92 1.13 1.12  -0.95   2.74 1.00     1932     1916
##    y1[21]   3.08   3.08 1.17 1.09   1.18   5.05 1.00     1890     1660
##    y1[22]  -1.36  -1.36 1.16 1.16  -3.20   0.59 1.00     1911     1916
##    y1[23]   0.20   0.22 1.14 1.15  -1.69   2.07 1.00     1596     1711
##    y1[24]  -2.56  -2.58 1.16 1.14  -4.52  -0.67 1.00     1621     1843
##    y1[25]  -1.20  -1.20 1.13 1.15  -3.04   0.66 1.00     2010     1780
##    y1[26]   2.83   2.86 1.22 1.14   0.77   4.79 1.00     1932     1834
##    y1[27]   0.28   0.27 1.13 1.11  -1.55   2.13 1.00     2092     1960
##    y1[28]   1.78   1.76 1.12 1.10  -0.08   3.64 1.00     2009     2012
##    y1[29]  -0.73  -0.71 1.15 1.12  -2.66   1.10 1.00     1811     1859
##    y1[30]  -2.34  -2.33 1.16 1.19  -4.25  -0.48 1.00     1546     1538
##    y1[31]   0.97   0.97 1.12 1.09  -0.89   2.78 1.00     2044     1917
##    y1[32]   0.08   0.10 1.14 1.13  -1.81   1.89 1.00     2023     1736
##    y1[33]   0.94   0.95 1.14 1.14  -1.00   2.79 1.00     1896     1840
##    y1[34]  -1.38  -1.37 1.15 1.17  -3.24   0.50 1.00     2059     1821
##    y1[35]   3.19   3.17 1.17 1.14   1.26   5.16 1.00     1986     1894
##    y1[36]   2.94   2.91 1.17 1.13   1.04   4.88 1.00     1918     1927
##    y1[37]  -1.18  -1.19 1.17 1.18  -3.08   0.78 1.00     2214     1970
##    y1[38]   0.46   0.47 1.13 1.12  -1.42   2.32 1.00     1873     2055
##    y1[39]   0.48   0.45 1.12 1.09  -1.34   2.37 1.00     1933     1926
##    y1[40]   0.50   0.48 1.15 1.09  -1.37   2.43 1.00     1696     2056
##    y1[41]   1.58   1.60 1.13 1.17  -0.31   3.34 1.00     1970     1968
##    y1[42]   2.78   2.78 1.13 1.14   0.93   4.59 1.00     1886     1776
##    y1[43]  -1.86  -1.86 1.20 1.19  -3.77   0.06 1.00     1417     1883
##    y1[44]   0.36   0.36 1.16 1.15  -1.51   2.33 1.00     1544     1885
##    y1[45]  -0.92  -0.91 1.14 1.14  -2.73   0.92 1.00     2105     1930
##    y1[46]   1.16   1.17 1.13 1.12  -0.70   3.00 1.00     2025     1959
##    y1[47]  -2.71  -2.74 1.20 1.16  -4.65  -0.66 1.00     1632     2008
##    y1[48]   2.91   2.91 1.23 1.22   0.86   4.96 1.00     1766     1793
##    y1[49]   3.86   3.87 1.18 1.15   1.86   5.70 1.00     1639     1841
##    y1[50]   1.24   1.27 1.18 1.14  -0.69   3.08 1.00     1673     1953
##    m1[1]   -1.39  -1.40 0.31 0.30  -1.89  -0.88 1.00     1698     1677
##    m1[2]   -1.53  -1.54 0.31 0.31  -2.04  -1.00 1.00     1724     1664
##    m1[3]    0.55   0.55 0.37 0.38  -0.06   1.14 1.00      879      907
##    m1[4]    1.11   1.11 0.28 0.28   0.63   1.56 1.00     1659     1369
##    m1[5]    1.52   1.52 0.29 0.28   1.05   1.98 1.00     1434     1504
##    m1[6]   -0.52  -0.52 0.27 0.27  -0.96  -0.06 1.00     1551     1594
##    m1[7]    1.01   1.01 0.28 0.28   0.54   1.46 1.00     1651     1295
##    m1[8]    1.32   1.32 0.36 0.37   0.73   1.92 1.00      862      891
##    m1[9]    3.19   3.20 0.34 0.33   2.63   3.73 1.00     2021     1591
##    m1[10]   2.57   2.58 0.31 0.30   2.06   3.07 1.00     1912     1619
##    m1[11]   3.47   3.48 0.36 0.34   2.88   4.03 1.00     2048     1671
##    m1[12]   3.40   3.40 0.39 0.38   2.75   4.01 1.00     1595     1614
##    m1[13]  -0.82  -0.83 0.28 0.28  -1.28  -0.36 1.00     1598     1587
##    m1[14]   2.22   2.22 0.42 0.42   1.51   2.91 1.00      666     1244
##    m1[15]   3.52   3.52 0.42 0.41   2.85   4.22 1.00     1059     1230
##    m1[16]  -1.51  -1.52 0.31 0.31  -2.02  -0.99 1.00     1721     1664
##    m1[17]  -0.69  -0.70 0.28 0.28  -1.15  -0.24 1.00     1578     1520
##    m1[18]   1.68   1.68 0.37 0.38   1.09   2.28 1.00      867      833
##    m1[19]   0.40   0.40 0.37 0.39  -0.22   1.00 1.00      887     1067
##    m1[20]   0.94   0.95 0.27 0.27   0.49   1.39 1.00     1423     1425
##    m1[21]   3.08   3.07 0.40 0.40   2.44   3.74 1.00     1007     1200
##    m1[22]  -1.33  -1.32 0.44 0.46  -2.04  -0.62 1.00     1060     1446
##    m1[23]   0.19   0.18 0.37 0.36  -0.40   0.77 1.00      540      828
##    m1[24]  -2.61  -2.61 0.43 0.43  -3.31  -1.89 1.00      697     1088
##    m1[25]  -1.17  -1.17 0.30 0.29  -1.65  -0.67 1.00     1660     1583
##    m1[26]   2.85   2.86 0.46 0.44   2.09   3.59 1.00      736     1332
##    m1[27]   0.29   0.28 0.29 0.29  -0.19   0.77 1.00     1577     1465
##    m1[28]   1.77   1.77 0.29 0.28   1.29   2.22 1.00     1756     1541
##    m1[29]  -0.73  -0.73 0.33 0.32  -1.27  -0.20 1.00     1593     1451
##    m1[30]  -2.36  -2.36 0.42 0.42  -3.05  -1.66 1.00      667      987
##    m1[31]   0.99   0.99 0.28 0.28   0.52   1.45 1.00     1650     1312
##    m1[32]   0.08   0.08 0.30 0.29  -0.41   0.56 1.00     1568     1420
##    m1[33]   0.96   0.96 0.27 0.27   0.51   1.41 1.00     1421     1425
##    m1[34]  -1.40  -1.40 0.36 0.35  -1.98  -0.81 1.00     1674     1595
##    m1[35]   3.22   3.23 0.34 0.33   2.66   3.76 1.00     2022     1617
##    m1[36]   2.95   2.96 0.33 0.31   2.42   3.48 1.00     1990     1493
##    m1[37]  -1.22  -1.22 0.35 0.34  -1.79  -0.66 1.00     1649     1581
##    m1[38]   0.47   0.47 0.26 0.26   0.03   0.91 1.00     1438     1542
##    m1[39]   0.50   0.49 0.29 0.28   0.02   0.96 1.00     1595     1444
##    m1[40]   0.47   0.47 0.26 0.26   0.03   0.91 1.00     1438     1542
##    m1[41]   1.59   1.60 0.29 0.29   1.11   2.06 1.00     1438     1564
##    m1[42]   2.76   2.75 0.39 0.39   2.14   3.40 1.00      957     1078
##    m1[43]  -1.89  -1.89 0.40 0.40  -2.54  -1.22 1.00      618     1010
##    m1[44]   0.36   0.35 0.37 0.37  -0.24   0.95 1.00      543      836
##    m1[45]  -0.89  -0.89 0.29 0.28  -1.35  -0.42 1.00     1610     1535
##    m1[46]   1.15   1.15 0.28 0.28   0.67   1.60 1.00     1663     1403
##    m1[47]  -2.74  -2.74 0.44 0.43  -3.45  -2.01 1.00      714     1101
##    m1[48]   2.93   2.94 0.46 0.45   2.16   3.67 1.00      745     1370
##    m1[49]   3.85   3.86 0.42 0.40   3.15   4.52 1.00     1631     1613
##    m1[50]   1.27   1.26 0.36 0.37   0.67   1.87 1.00      861      891

tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,xa,y,col=ca),
             qplot(data=tb,xa,y,col=cb),
             qplot(data=tb,xb,y,col=ca),
             qplot(data=tb,xb,y,col=cb))

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)

fn(f0)
## Initial log joint probability = -847.296 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -48.2121   0.000112663     0.0012299           1           1       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -50.93 -50.60 1.81 1.70 -54.47 -48.60 1.00      870     1288
##    b[1]    -0.02  -0.01 0.41 0.38  -0.70   0.65 1.00      660      812
##    b[2]     1.25   1.25 0.20 0.20   0.92   1.57 1.00     2716     1608
##    b[3]     1.09   1.10 0.23 0.24   0.72   1.46 1.00     2081     1523
##    b[4]     2.04   2.03 0.48 0.46   1.23   2.83 1.00      733      785
##    b[5]     1.50   1.48 0.51 0.50   0.68   2.33 1.01      701      856
##    s        1.72   1.71 0.19 0.17   1.44   2.04 1.00     1965     1437
##    y1[1]   -3.26  -3.22 1.83 1.85  -6.20  -0.30 1.00     1881     1834
##    y1[2]    0.65   0.66 1.80 1.72  -2.22   3.72 1.00     2032     1895
##    y1[3]    1.79   1.81 1.79 1.73  -1.20   4.69 1.00     1963     1641
##    y1[4]   -0.81  -0.83 1.87 1.80  -3.94   2.26 1.00     1885     1920
##    y1[5]    3.20   3.16 1.80 1.75   0.25   6.21 1.00     2000     1969
##    y1[6]    4.87   4.86 1.74 1.67   2.03   7.79 1.00     1963     1996
##    y1[7]   -2.14  -2.07 1.80 1.81  -5.13   0.82 1.00     1784     1932
##    y1[8]    0.97   0.97 1.81 1.76  -1.89   4.00 1.00     1931     1892
##    y1[9]    3.02   2.95 1.82 1.79   0.06   5.96 1.00     1899     1970
##    y1[10]   0.41   0.40 1.84 1.82  -2.67   3.49 1.00     1959     1954
##    y1[11]  -2.39  -2.43 1.83 1.79  -5.43   0.67 1.00     1957     1953
##    y1[12]   3.71   3.77 1.87 1.83   0.72   6.67 1.00     2056     2019
##    y1[13]   5.19   5.20 1.87 1.83   2.11   8.31 1.00     1965     1971
##    y1[14]   1.13   1.12 1.74 1.70  -1.64   4.07 1.00     1966     1819
##    y1[15]  -0.59  -0.60 1.84 1.79  -3.60   2.58 1.00     2039     1812
##    y1[16]  -2.52  -2.53 1.84 1.75  -5.49   0.52 1.00     1761     1766
##    y1[17]   4.12   4.13 1.80 1.72   1.18   7.08 1.00     1521     1629
##    y1[18]   0.58   0.60 1.79 1.78  -2.39   3.49 1.00     2013     1972
##    y1[19]   1.29   1.31 1.80 1.82  -1.74   4.19 1.00     1944     2015
##    y1[20]  -0.18  -0.18 1.78 1.79  -3.07   2.86 1.00     1971     1951
##    y1[21]  -3.08  -3.05 1.76 1.70  -5.93  -0.25 1.00     1840     1573
##    y1[22]   2.36   2.35 1.91 1.90  -0.80   5.48 1.00     1894     2014
##    y1[23]   5.69   5.67 1.88 1.85   2.59   8.73 1.00     1946     1863
##    y1[24]   3.12   3.12 1.84 1.81   0.11   6.14 1.00     1704     1963
##    y1[25]  -0.08  -0.06 1.77 1.69  -3.04   2.81 1.00     1738     1878
##    y1[26]  -0.45  -0.44 1.79 1.75  -3.35   2.50 1.00     1827     1717
##    y1[27]  -0.12  -0.08 1.82 1.77  -3.10   2.82 1.00     1507     1890
##    y1[28]   2.71   2.72 1.80 1.76  -0.26   5.69 1.00     1764     1664
##    y1[29]   1.89   1.88 1.81 1.81  -1.18   4.80 1.00     1778     1709
##    y1[30]   5.94   5.97 1.89 1.83   2.79   9.04 1.00     2017     1966
##    y1[31]   2.21   2.19 1.89 1.79  -0.85   5.35 1.00     1904     2047
##    y1[32]   3.99   4.00 1.87 1.91   0.96   7.04 1.00     2076     1974
##    y1[33]  -1.19  -1.20 1.76 1.75  -4.05   1.74 1.00     2158     2093
##    y1[34]  -1.14  -1.14 1.88 1.87  -4.12   1.93 1.00     2050     1913
##    y1[35]   2.30   2.31 1.82 1.83  -0.67   5.26 1.00     1953     1879
##    y1[36]   5.96   5.95 1.82 1.79   2.96   8.97 1.00     2173     1880
##    y1[37]   2.09   2.11 1.77 1.76  -0.89   4.97 1.00     1989     1944
##    y1[38]   5.35   5.34 1.77 1.80   2.49   8.26 1.00     1934     1919
##    y1[39]   0.95   0.92 1.89 1.76  -2.08   4.11 1.00     1997     2014
##    y1[40]  -0.48  -0.47 1.76 1.81  -3.36   2.40 1.00     1602     1733
##    y1[41]   0.87   0.81 1.80 1.76  -2.14   3.82 1.00     1582     1891
##    y1[42]  -0.26  -0.27 1.78 1.72  -3.11   2.67 1.00     2019     1822
##    y1[43]  -2.14  -2.11 1.78 1.73  -5.22   0.76 1.00     1864     1833
##    y1[44]   5.06   5.05 1.77 1.75   2.17   7.93 1.00     1920     1991
##    y1[45]   0.30   0.31 1.80 1.74  -2.58   3.34 1.00     2116     1918
##    y1[46]  -3.69  -3.72 1.82 1.76  -6.73  -0.70 1.00     2000     1687
##    y1[47]   2.41   2.43 1.79 1.76  -0.45   5.44 1.00     1526     1832
##    y1[48]   0.12   0.07 1.81 1.79  -2.83   3.05 1.00     1762     1787
##    y1[49]   4.89   4.91 1.78 1.72   1.94   7.69 1.00     1937     1934
##    y1[50]   1.49   1.51 1.87 1.74  -1.64   4.67 1.00     1883     1722
##    m1[1]   -3.21  -3.22 0.51 0.52  -4.04  -2.33 1.00     1301     1239
##    m1[2]    0.68   0.68 0.61 0.60  -0.32   1.72 1.00     1366     1273
##    m1[3]    1.84   1.85 0.49 0.49   1.07   2.63 1.00     1385     1536
##    m1[4]   -0.77  -0.79 0.57 0.57  -1.70   0.19 1.00     1405     1431
##    m1[5]    3.29   3.29 0.49 0.48   2.49   4.07 1.00     1267     1419
##    m1[6]    4.92   4.91 0.47 0.45   4.17   5.68 1.00     1248     1315
##    m1[7]   -2.11  -2.11 0.51 0.50  -2.96  -1.26 1.00     1219     1311
##    m1[8]    0.92   0.92 0.49 0.49   0.13   1.74 1.00     1167     1335
##    m1[9]    2.97   2.97 0.60 0.59   1.95   3.96 1.00      956     1215
##    m1[10]   0.41   0.41 0.67 0.65  -0.71   1.49 1.00     1740     1474
##    m1[11]  -2.39  -2.39 0.53 0.52  -3.27  -1.54 1.00     1181     1391
##    m1[12]   3.73   3.72 0.55 0.56   2.84   4.61 1.00     1294     1290
##    m1[13]   5.23   5.22 0.53 0.52   4.37   6.10 1.00     1433     1325
##    m1[14]   1.09   1.08 0.46 0.46   0.37   1.83 1.00     1164     1267
##    m1[15]  -0.56  -0.57 0.62 0.61  -1.59   0.46 1.00     1565     1520
##    m1[16]  -2.55  -2.55 0.50 0.49  -3.37  -1.71 1.00     1134     1297
##    m1[17]   4.14   4.14 0.56 0.54   3.21   5.06 1.00     1469     1433
##    m1[18]   0.61   0.60 0.44 0.42  -0.11   1.34 1.00     1257     1307
##    m1[19]   1.29   1.28 0.43 0.43   0.59   2.00 1.00     1188     1245
##    m1[20]  -0.22  -0.23 0.47 0.45  -0.98   0.55 1.00     1419     1343
##    m1[21]  -3.03  -3.03 0.52 0.51  -3.84  -2.15 1.00     1310     1305
##    m1[22]   2.39   2.39 0.69 0.68   1.22   3.51 1.00     1352     1332
##    m1[23]   5.77   5.78 0.66 0.66   4.67   6.85 1.00     1759     1655
##    m1[24]   3.12   3.14 0.64 0.63   2.04   4.16 1.00      938     1287
##    m1[25]   0.00   0.00 0.50 0.49  -0.83   0.82 1.00      956     1007
##    m1[26]  -0.45  -0.47 0.47 0.47  -1.22   0.32 1.00     1417     1340
##    m1[27]  -0.11  -0.11 0.66 0.67  -1.17   0.98 1.00     1273     1537
##    m1[28]   2.67   2.66 0.49 0.49   1.87   3.48 1.00     1108     1402
##    m1[29]   1.86   1.86 0.48 0.48   1.09   2.65 1.00     1171     1267
##    m1[30]   5.89   5.90 0.71 0.73   4.74   7.04 1.00     1374     1525
##    m1[31]   2.21   2.21 0.65 0.66   1.16   3.24 1.00     1875     1636
##    m1[32]   4.04   4.04 0.51 0.52   3.20   4.87 1.00     1242     1242
##    m1[33]  -1.18  -1.20 0.55 0.53  -2.09  -0.28 1.00     1784     1410
##    m1[34]  -1.13  -1.15 0.55 0.53  -2.01  -0.23 1.00     1554     1305
##    m1[35]   2.35   2.34 0.58 0.56   1.38   3.30 1.00     1202     1183
##    m1[36]   5.97   5.95 0.52 0.50   5.12   6.85 1.00     1524     1378
##    m1[37]   2.10   2.10 0.45 0.45   1.39   2.83 1.00     1119     1249
##    m1[38]   5.30   5.30 0.49 0.47   4.51   6.11 1.00     1338     1343
##    m1[39]   0.91   0.90 0.62 0.61  -0.15   1.90 1.00     1649     1431
##    m1[40]  -0.44  -0.44 0.43 0.40  -1.14   0.25 1.00      708      864
##    m1[41]   0.87   0.88 0.58 0.56  -0.08   1.83 1.00     1101     1327
##    m1[42]  -0.31  -0.32 0.46 0.46  -1.09   0.46 1.00      897      977
##    m1[43]  -2.05  -2.06 0.44 0.42  -2.78  -1.33 1.00      914     1116
##    m1[44]   5.08   5.06 0.48 0.47   4.30   5.89 1.00     1319     1350
##    m1[45]   0.27   0.27 0.44 0.44  -0.46   1.00 1.00     1275     1301
##    m1[46]  -3.66  -3.67 0.55 0.53  -4.55  -2.74 1.00     1460     1278
##    m1[47]   2.38   2.38 0.59 0.57   1.38   3.33 1.00      874     1222
##    m1[48]   0.11   0.11 0.65 0.64  -0.95   1.21 1.00     1294     1563
##    m1[49]   4.84   4.82 0.49 0.49   4.04   5.68 1.00     1314     1327
##    m1[50]   1.51   1.50 0.72 0.71   0.31   2.73 1.00     1447     1380

fn(f1)
## Initial log joint probability = -116.493 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       27      -45.1234   0.000147958   0.000662463           1           1       31    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -48.59 -48.23 2.09 1.84 -52.59 -45.95 1.01      695      957
##    b[1]    -0.51  -0.52 0.46 0.48  -1.24   0.21 1.01      655      871
##    b[2]     1.29   1.29 0.19 0.19   0.98   1.59 1.00     2970     1618
##    b[3]     1.04   1.04 0.22 0.22   0.66   1.39 1.00     2333     1506
##    b[4]     3.02   3.02 0.63 0.63   2.02   4.04 1.01      570     1035
##    b[5]     2.58   2.57 0.69 0.71   1.50   3.71 1.01      571      761
##    b[6]    -2.19  -2.16 0.97 0.98  -3.81  -0.60 1.02      495      590
##    s        1.65   1.64 0.19 0.18   1.38   1.97 1.00     1559     1317
##    y1[1]   -3.68  -3.66 1.71 1.71  -6.54  -0.84 1.00     1791     1770
##    y1[2]    1.40   1.34 1.76 1.74  -1.41   4.33 1.00     1647     1932
##    y1[3]    2.45   2.45 1.77 1.76  -0.42   5.32 1.00     2145     2055
##    y1[4]   -0.24  -0.23 1.73 1.68  -3.12   2.63 1.00     1890     1925
##    y1[5]    3.91   3.89 1.79 1.78   0.94   6.89 1.00     1811     1973
##    y1[6]    4.31   4.34 1.75 1.69   1.43   7.18 1.00     1895     2009
##    y1[7]   -2.46  -2.45 1.75 1.77  -5.39   0.34 1.00     1896     1746
##    y1[8]    1.58   1.57 1.75 1.72  -1.27   4.38 1.00     1853     1919
##    y1[9]    2.52   2.57 1.76 1.80  -0.36   5.38 1.00     1805     1909
##    y1[10]   0.88   0.83 1.80 1.73  -1.99   3.91 1.00     1875     1916
##    y1[11]  -2.94  -2.92 1.74 1.75  -5.72  -0.07 1.00     1906     1861
##    y1[12]   3.22   3.23 1.76 1.75   0.37   6.18 1.00     2011     2083
##    y1[13]   4.67   4.65 1.79 1.78   1.81   7.63 1.00     1988     1945
##    y1[14]   1.67   1.66 1.73 1.69  -1.21   4.52 1.00     1937     1895
##    y1[15]  -0.05  -0.07 1.75 1.71  -2.95   2.79 1.00     2034     2011
##    y1[16]  -3.09  -3.13 1.75 1.70  -5.87  -0.27 1.00     1724     1811
##    y1[17]   4.84   4.81 1.77 1.66   1.98   7.83 1.00     1980     1971
##    y1[18]   1.09   1.15 1.73 1.70  -1.80   3.84 1.00     2124     1872
##    y1[19]   1.79   1.80 1.67 1.67  -0.96   4.47 1.00     2175     2025
##    y1[20]   0.34   0.32 1.69 1.66  -2.41   3.07 1.00     1635     1706
##    y1[21]  -3.48  -3.45 1.71 1.68  -6.33  -0.71 1.00     1830     1922
##    y1[22]   2.66   2.65 1.79 1.74  -0.35   5.68 1.00     1918     1929
##    y1[23]   6.31   6.35 1.82 1.74   3.29   9.33 1.00     2054     1969
##    y1[24]   2.59   2.63 1.81 1.80  -0.45   5.47 1.00     1743     1995
##    y1[25]  -0.39  -0.40 1.75 1.75  -3.19   2.53 1.00     2099     1892
##    y1[26]   0.04   0.07 1.75 1.64  -2.86   2.90 1.00     1886     1827
##    y1[27]  -0.67  -0.65 1.81 1.75  -3.64   2.31 1.00     1946     1972
##    y1[28]   2.06   2.07 1.73 1.66  -0.76   4.90 1.00     1939     1821
##    y1[29]   2.45   2.46 1.75 1.74  -0.43   5.37 1.00     1512     1636
##    y1[30]   6.35   6.31 1.79 1.72   3.40   9.21 1.00     1990     1807
##    y1[31]   2.92   2.90 1.84 1.75  -0.04   6.05 1.00     1918     2059
##    y1[32]   3.53   3.48 1.70 1.71   0.85   6.28 1.00     1697     1892
##    y1[33]  -0.58  -0.53 1.77 1.76  -3.46   2.36 1.00     2018     2015
##    y1[34]  -0.68  -0.67 1.76 1.77  -3.65   2.21 1.00     2023     1796
##    y1[35]   2.74   2.78 1.74 1.71  -0.16   5.63 1.00     1915     1779
##    y1[36]   5.38   5.40 1.72 1.69   2.50   8.24 1.00     1954     1845
##    y1[37]   2.69   2.66 1.74 1.74  -0.29   5.46 1.00     1835     1762
##    y1[38]   4.66   4.67 1.75 1.76   1.80   7.50 1.00     2082     2007
##    y1[39]   1.40   1.40 1.81 1.75  -1.59   4.29 1.00     1985     1894
##    y1[40]  -0.96  -0.99 1.75 1.79  -3.78   1.78 1.00     1815     1735
##    y1[41]   0.59   0.60 1.74 1.69  -2.26   3.50 1.00     1916     1971
##    y1[42]  -0.75  -0.73 1.70 1.62  -3.48   1.99 1.00     1927     1980
##    y1[43]  -2.52  -2.50 1.72 1.70  -5.38   0.25 1.00     1884     2008
##    y1[44]   4.40   4.41 1.73 1.75   1.53   7.25 1.00     2002     1960
##    y1[45]   0.78   0.79 1.70 1.70  -1.91   3.60 1.00     1776     1724
##    y1[46]  -4.27  -4.28 1.75 1.70  -7.08  -1.39 1.00     1906     1824
##    y1[47]   1.84   1.85 1.74 1.63  -1.04   4.70 1.00     1848     1819
##    y1[48]  -0.53  -0.52 1.77 1.73  -3.46   2.39 1.00     1948     1932
##    y1[49]   4.12   4.16 1.74 1.77   1.25   6.90 1.00     2070     1786
##    y1[50]   1.85   1.80 1.77 1.70  -0.97   4.78 1.00     1872     1827
##    m1[1]   -3.69  -3.69 0.52 0.52  -4.52  -2.82 1.01      960     1310
##    m1[2]    1.39   1.37 0.69 0.69   0.33   2.57 1.01     1134     1222
##    m1[3]    2.43   2.43 0.54 0.53   1.55   3.33 1.00     1431     1533
##    m1[4]   -0.24  -0.25 0.59 0.59  -1.17   0.73 1.00     1305     1350
##    m1[5]    3.91   3.90 0.57 0.57   3.00   4.88 1.00     1135     1235
##    m1[6]    4.30   4.29 0.54 0.52   3.42   5.16 1.00     1578     1204
##    m1[7]   -2.51  -2.52 0.53 0.54  -3.38  -1.65 1.00      996     1622
##    m1[8]    1.57   1.55 0.57 0.57   0.68   2.50 1.01     1024     1149
##    m1[9]    2.52   2.52 0.63 0.64   1.48   3.52 1.01     1075     1063
##    m1[10]   0.84   0.82 0.67 0.65  -0.20   1.93 1.00     1721     1487
##    m1[11]  -2.96  -2.95 0.54 0.53  -3.85  -2.07 1.01      854     1000
##    m1[12]   3.22   3.21 0.60 0.60   2.26   4.20 1.00     1834     1507
##    m1[13]   4.69   4.69 0.58 0.57   3.74   5.67 1.00     1820     1468
##    m1[14]   1.65   1.64 0.51 0.51   0.83   2.49 1.00     1090     1419
##    m1[15]  -0.08  -0.09 0.62 0.61  -1.04   0.94 1.00     1483     1443
##    m1[16]  -3.09  -3.09 0.51 0.50  -3.92  -2.26 1.01      823      969
##    m1[17]   4.79   4.79 0.64 0.62   3.77   5.88 1.00     1262     1342
##    m1[18]   1.15   1.15 0.48 0.47   0.37   1.95 1.00     1351     1590
##    m1[19]   1.83   1.84 0.48 0.46   1.06   2.62 1.00     1322     1433
##    m1[20]   0.31   0.32 0.50 0.50  -0.51   1.16 1.00     1461     1520
##    m1[21]  -3.48  -3.48 0.53 0.54  -4.33  -2.59 1.00      990     1515
##    m1[22]   2.74   2.74 0.67 0.64   1.65   3.81 1.00     1667     1665
##    m1[23]   6.34   6.33 0.71 0.68   5.16   7.49 1.00     1606     1455
##    m1[24]   2.60   2.60 0.67 0.68   1.50   3.67 1.01     1034     1055
##    m1[25]  -0.39  -0.39 0.52 0.53  -1.25   0.47 1.00      937     1344
##    m1[26]   0.06   0.06 0.50 0.49  -0.74   0.89 1.00     1504     1572
##    m1[27]  -0.70  -0.69 0.69 0.66  -1.88   0.41 1.00     1809     1540
##    m1[28]   2.01   2.01 0.55 0.56   1.09   2.90 1.00     1531     1405
##    m1[29]   2.51   2.49 0.57 0.55   1.64   3.43 1.01     1024     1110
##    m1[30]   6.39   6.40 0.71 0.69   5.20   7.48 1.00     1612     1519
##    m1[31]   2.88   2.87 0.70 0.68   1.73   4.03 1.00     1710     1796
##    m1[32]   3.51   3.50 0.57 0.57   2.58   4.44 1.00     1747     1638
##    m1[33]  -0.62  -0.61 0.58 0.58  -1.59   0.36 1.00     1632     1417
##    m1[34]  -0.68  -0.70 0.54 0.52  -1.54   0.25 1.00     1715     1509
##    m1[35]   2.75   2.74 0.57 0.55   1.81   3.68 1.00     1492     1553
##    m1[36]   5.37   5.38 0.58 0.55   4.44   6.34 1.00     1735     1369
##    m1[37]   2.71   2.71 0.53 0.53   1.88   3.59 1.00     1022     1179
##    m1[38]   4.71   4.70 0.55 0.53   3.81   5.63 1.00     1653     1398
##    m1[39]   1.36   1.34 0.63 0.61   0.37   2.42 1.00     1597     1360
##    m1[40]  -0.96  -0.97 0.46 0.47  -1.72  -0.22 1.01      657      975
##    m1[41]   0.51   0.52 0.60 0.61  -0.48   1.49 1.00     1190     1713
##    m1[42]  -0.72  -0.73 0.49 0.51  -1.54   0.08 1.00      853     1280
##    m1[43]  -2.53  -2.53 0.46 0.47  -3.27  -1.76 1.01      745     1093
##    m1[44]   4.44   4.43 0.55 0.53   3.54   5.34 1.00     1583     1288
##    m1[45]   0.80   0.81 0.48 0.46   0.02   1.59 1.00     1386     1586
##    m1[46]  -4.16  -4.16 0.55 0.55  -5.04  -3.26 1.01     1024     1226
##    m1[47]   1.84   1.85 0.63 0.65   0.82   2.85 1.01      913      918
##    m1[48]  -0.54  -0.53 0.68 0.67  -1.68   0.55 1.00     1781     1536
##    m1[49]   4.18   4.17 0.56 0.54   3.26   5.09 1.00     1569     1325
##    m1[50]   1.85   1.84 0.69 0.64   0.73   2.97 1.00     1777     1671

fn(f2)
## Initial log joint probability = -1164.03 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       28      -18.0098   7.85082e-05   0.000763031           1           1       35    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -22.53 -22.19 2.16 2.00 -26.56 -19.71 1.00      745     1205
##    b[1]     0.01   0.00 0.28 0.28  -0.44   0.48 1.00      673      985
##    b[2]     1.24   1.24 0.11 0.11   1.05   1.42 1.00     2311     1595
##    b[3]     1.04   1.04 0.13 0.13   0.84   1.26 1.00     2457     1955
##    b[4]     2.37   2.38 0.38 0.38   1.74   2.99 1.01      524      723
##    b[5]     1.86   1.86 0.40 0.40   1.19   2.52 1.00      450      757
##    b[6]    -0.91  -0.91 0.10 0.10  -1.07  -0.74 1.00     2037     1465
##    b[7]    -1.03  -1.05 0.58 0.59  -1.93  -0.04 1.00      398      620
##    s        0.97   0.96 0.11 0.11   0.81   1.17 1.00     1958     1703
##    y1[1]   -4.81  -4.80 1.01 0.97  -6.53  -3.18 1.00     2124     2002
##    y1[2]    2.39   2.42 1.03 1.02   0.65   4.08 1.00     1908     1932
##    y1[3]    3.16   3.17 1.05 0.99   1.37   4.88 1.00     1772     1627
##    y1[4]   -0.74  -0.73 1.04 1.08  -2.43   0.98 1.00     1970     1583
##    y1[5]    3.40   3.43 1.04 1.02   1.75   5.05 1.00     1872     1792
##    y1[6]    4.25   4.28 1.04 1.02   2.54   5.98 1.00     1979     2035
##    y1[7]   -2.20  -2.22 1.00 1.02  -3.78  -0.56 1.00     1891     2009
##    y1[8]    1.62   1.62 1.04 1.04  -0.07   3.35 1.00     1574     1611
##    y1[9]    1.80   1.79 1.05 0.99   0.10   3.57 1.00     1779     1732
##    y1[10]   2.85   2.84 1.06 1.01   1.10   4.65 1.00     2110     1930
##    y1[11]  -2.49  -2.48 1.04 1.07  -4.20  -0.83 1.00     1645     1854
##    y1[12]   4.69   4.67 1.05 1.05   3.00   6.43 1.00     2053     1820
##    y1[13]   5.34   5.34 1.01 1.01   3.67   6.98 1.00     1879     1884
##    y1[14]   1.58   1.56 1.03 1.00  -0.15   3.24 1.00     1781     1702
##    y1[15]   0.39   0.41 1.01 0.98  -1.27   2.00 1.00     1871     1860
##    y1[16]  -3.14  -3.12 1.02 0.99  -4.88  -1.52 1.00     1450     1868
##    y1[17]   4.15   4.13 1.04 1.00   2.44   5.92 1.00     2023     1759
##    y1[18]   0.89   0.89 1.02 0.97  -0.77   2.60 1.00     1805     1777
##    y1[19]   1.84   1.85 1.01 0.99   0.17   3.51 1.00     1883     1873
##    y1[20]  -0.53  -0.51 1.05 1.03  -2.26   1.17 1.00     1585     1487
##    y1[21]  -4.34  -4.35 1.02 0.99  -5.97  -2.62 1.00     1857     1876
##    y1[22]   4.47   4.46 1.06 1.03   2.79   6.22 1.00     2081     1972
##    y1[23]   2.98   3.00 1.14 1.16   1.06   4.83 1.00     1584     1643
##    y1[24]   1.40   1.41 1.06 1.05  -0.33   3.12 1.00     2020     2016
##    y1[25]   1.01   0.98 1.04 1.03  -0.66   2.79 1.00     1715     1867
##    y1[26]  -1.02  -1.01 1.03 1.02  -2.69   0.65 1.00     1853     1621
##    y1[27]  -2.49  -2.46 1.09 1.07  -4.28  -0.72 1.00     2038     2010
##    y1[28]   2.42   2.38 1.04 0.99   0.74   4.19 1.00     2000     2016
##    y1[29]   2.72   2.71 1.02 1.00   1.04   4.43 1.00     1877     1880
##    y1[30]   3.78   3.79 1.09 1.06   1.96   5.52 1.00     1761     1787
##    y1[31]   5.63   5.63 1.08 1.07   3.87   7.44 1.00     1839     1854
##    y1[32]   4.57   4.59 1.02 0.98   2.89   6.21 1.00     1938     1798
##    y1[33]  -2.24  -2.26 1.07 1.05  -4.04  -0.50 1.00     1748     1947
##    y1[34]  -2.00  -1.99 1.01 1.02  -3.63  -0.34 1.00     1719     1729
##    y1[35]   3.42   3.43 1.03 0.98   1.75   5.10 1.00     1901     1965
##    y1[36]   4.78   4.81 1.03 1.05   3.05   6.47 1.00     1904     1890
##    y1[37]   2.52   2.51 1.02 1.01   0.86   4.23 1.00     1639     1916
##    y1[38]   4.61   4.61 1.03 0.99   2.97   6.26 1.00     1897     1972
##    y1[39]   2.91   2.92 1.03 1.02   1.22   4.63 1.00     2131     1907
##    y1[40]  -0.36  -0.35 0.97 0.96  -1.93   1.24 1.00     1533     1603
##    y1[41]   2.60   2.61 1.05 1.07   0.89   4.34 1.00     1719     1936
##    y1[42]   0.33   0.33 1.02 1.01  -1.33   1.96 1.00     1833     1853
##    y1[43]  -2.64  -2.63 0.99 0.96  -4.31  -1.07 1.00     2122     1650
##    y1[44]   4.31   4.30 1.04 1.01   2.64   6.02 1.00     1864     1996
##    y1[45]   0.28   0.28 1.02 1.00  -1.35   1.98 1.00     1836     1839
##    y1[46]  -5.64  -5.66 1.06 1.04  -7.37  -3.91 1.00     1976     2054
##    y1[47]   1.42   1.40 1.03 1.04  -0.20   3.15 1.00     1385     1805
##    y1[48]  -1.76  -1.74 1.06 1.01  -3.46  -0.04 1.00     1781     1834
##    y1[49]   4.29   4.29 1.02 1.01   2.68   5.99 1.00     1884     1773
##    y1[50]   4.29   4.32 1.11 1.12   2.50   6.07 1.00     1951     1896
##    m1[1]   -4.79  -4.80 0.33 0.32  -5.33  -4.24 1.00     1428     1627
##    m1[2]    2.37   2.37 0.39 0.40   1.73   3.03 1.00      966     1626
##    m1[3]    3.15   3.15 0.32 0.32   2.64   3.67 1.00     1091     1566
##    m1[4]   -0.73  -0.72 0.34 0.33  -1.29  -0.17 1.00      991     1296
##    m1[5]    3.40   3.41 0.32 0.31   2.83   3.92 1.00      721     1073
##    m1[6]    4.28   4.29 0.31 0.30   3.78   4.79 1.00     1428     1579
##    m1[7]   -2.19  -2.20 0.31 0.30  -2.69  -1.67 1.00     1161     1240
##    m1[8]    1.61   1.61 0.31 0.30   1.10   2.13 1.00      707     1300
##    m1[9]    1.82   1.82 0.39 0.37   1.18   2.43 1.00     1080     1164
##    m1[10]   2.85   2.85 0.45 0.45   2.11   3.59 1.00     2124     1714
##    m1[11]  -2.51  -2.52 0.33 0.33  -3.04  -1.97 1.00      810     1101
##    m1[12]   4.74   4.75 0.36 0.36   4.14   5.34 1.00     1624     1569
##    m1[13]   5.34   5.35 0.32 0.32   4.80   5.88 1.00     1632     1584
##    m1[14]   1.54   1.54 0.28 0.28   1.06   1.99 1.00      786     1324
##    m1[15]   0.40   0.40 0.36 0.36  -0.19   1.00 1.00     1417     1427
##    m1[16]  -3.13  -3.14 0.31 0.30  -3.62  -2.61 1.00      856     1284
##    m1[17]   4.11   4.12 0.37 0.34   3.47   4.70 1.00      800     1205
##    m1[18]   0.84   0.84 0.29 0.28   0.35   1.31 1.01      887     1386
##    m1[19]   1.85   1.85 0.28 0.27   1.38   2.32 1.01      884     1343
##    m1[20]  -0.55  -0.54 0.32 0.31  -1.08  -0.04 1.01      949     1418
##    m1[21]  -4.36  -4.37 0.32 0.32  -4.88  -3.82 1.00     1440     1543
##    m1[22]   4.48   4.48 0.45 0.45   3.76   5.22 1.00     1937     1585
##    m1[23]   2.99   3.01 0.56 0.55   2.05   3.88 1.00     1145     1304
##    m1[24]   1.42   1.43 0.42 0.41   0.73   2.10 1.00     1103     1272
##    m1[25]   1.03   1.03 0.34 0.34   0.48   1.59 1.00     1008     1297
##    m1[26]  -1.03  -1.03 0.33 0.32  -1.57  -0.50 1.01      991     1507
##    m1[27]  -2.50  -2.50 0.46 0.45  -3.27  -1.75 1.00     1903     1451
##    m1[28]   2.41   2.42 0.33 0.31   1.85   2.93 1.00     1269     1400
##    m1[29]   2.70   2.70 0.31 0.30   2.18   3.22 1.00      726     1285
##    m1[30]   3.78   3.79 0.52 0.50   2.90   4.63 1.01     1186     1448
##    m1[31]   5.60   5.60 0.49 0.50   4.82   6.40 1.00     1970     1856
##    m1[32]   4.58   4.59 0.33 0.32   4.04   5.11 1.00     1541     1491
##    m1[33]  -2.20  -2.21 0.39 0.38  -2.84  -1.55 1.00     1096     1495
##    m1[34]  -1.99  -1.99 0.36 0.37  -2.57  -1.42 1.01     1221     1400
##    m1[35]   3.43   3.42 0.35 0.36   2.86   4.01 1.01     1476     1302
##    m1[36]   4.79   4.79 0.34 0.34   4.25   5.36 1.00     1669     1558
##    m1[37]   2.53   2.54 0.29 0.28   2.03   3.00 1.00      678     1053
##    m1[38]   4.61   4.60 0.31 0.31   4.10   5.12 1.00     1527     1506
##    m1[39]   2.92   2.92 0.40 0.40   2.27   3.57 1.00     1942     1658
##    m1[40]  -0.35  -0.36 0.29 0.28  -0.82   0.12 1.00      663      941
##    m1[41]   2.61   2.61 0.42 0.42   1.94   3.29 1.00     1202     1428
##    m1[42]   0.38   0.37 0.31 0.32  -0.12   0.89 1.00      922     1192
##    m1[43]  -2.65  -2.66 0.28 0.28  -3.09  -2.19 1.00      846     1327
##    m1[44]   4.32   4.33 0.32 0.31   3.81   4.85 1.00     1443     1563
##    m1[45]   0.25   0.25 0.29 0.29  -0.25   0.72 1.01      901     1420
##    m1[46]  -5.64  -5.65 0.36 0.36  -6.23  -5.03 1.00     1631     1545
##    m1[47]   1.45   1.45 0.37 0.36   0.85   2.06 1.00      920     1206
##    m1[48]  -1.76  -1.76 0.43 0.42  -2.48  -1.04 1.00     1770     1404
##    m1[49]   4.26   4.26 0.33 0.32   3.73   4.79 1.00     1401     1480
##    m1[50]   4.31   4.30 0.51 0.52   3.48   5.13 1.00     2083     1730

generalized linear regression

log normal regression

objective variable [0,Infinity)

# of samples is N,  
log mi=sum(bj*xji),j=[0,K],i=[1,N]  
log yi~N(mi,s)  

ex5-2.stan

log normal regression

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> s;
}
model{
  vector[N] m=X*b;
  y~lognormal(m,s);
}
generated quantities{
  vector[N] y1;
  vector[N] m1=X*b;
  for(i in 1:N){
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
          y=rnorm(n,log(x1+x2),1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex5-2.stan')

mle=mdl$optimize(data=data)  # it sometimes occur error and stop process
## Initial log joint probability = -1103.3 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       32      -11.1287   4.18458e-05     0.0013495      0.9083      0.9083       46    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -11.13
##    b[1]       0.81
##    b[2]       0.08
##    b[3]      -0.07
##    s          0.42
##    y1[1]      1.88
##    y1[2]      6.51
##    y1[3]      5.03
##    y1[4]      3.10
##    y1[5]      3.86
##    y1[6]      2.22
##    y1[7]      3.00
##    y1[8]      4.03
##    y1[9]      2.58
##    y1[10]     2.96
##    y1[11]     6.59
##    y1[12]     1.35
##    y1[13]     2.04
##    y1[14]     2.06
##    y1[15]     2.73
##    y1[16]     3.81
##    y1[17]     3.13
##    y1[18]     1.81
##    y1[19]     0.80
##    y1[20]     6.89
##    m1[1]      0.87
##    m1[2]      1.35
##    m1[3]      1.01
##    m1[4]      0.72
##    m1[5]      0.90
##    m1[6]      0.72
##    m1[7]      1.00
##    m1[8]      1.15
##    m1[9]      0.88
##    m1[10]     0.78
##    m1[11]     1.31
##    m1[12]     0.38
##    m1[13]     1.28
##    m1[14]     0.58
##    m1[15]     0.59
##    m1[16]     1.11
##    m1[17]     1.19
##    m1[18]     0.41
##    m1[19]     0.21
##    m1[20]     1.33
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.30 -13.86 1.74 1.39 -17.72 -12.40 1.01      386      567
##    b[1]     0.81   0.81 0.34 0.34   0.23   1.38 1.02      283      360
##    b[2]     0.08   0.08 0.04 0.04   0.01   0.15 1.01      466      516
##    b[3]    -0.07  -0.07 0.04 0.04  -0.14   0.00 1.01      647      614
##    s        0.50   0.49 0.10 0.09   0.37   0.71 1.01      486      648
##    y1[1]    2.73   2.38 1.68 1.21   0.94   5.69 1.00     2048     1881
##    y1[2]    4.54   3.86 2.82 1.98   1.59   9.75 1.00     1774     1890
##    y1[3]    3.12   2.74 1.79 1.38   1.12   6.33 1.00     1756     1618
##    y1[4]    2.38   2.06 1.61 1.03   0.83   4.98 1.00     2039     1864
##    y1[5]    2.83   2.41 1.80 1.20   0.99   5.99 1.00     1883     1935
##    y1[6]    2.38   2.10 1.41 1.06   0.82   4.94 1.00     1484     1868
##    y1[7]    3.10   2.71 1.79 1.33   1.15   6.17 1.00     1860     1414
##    y1[8]    3.60   3.11 2.19 1.56   1.28   7.42 1.00     1984     1666
##    y1[9]    2.87   2.45 1.87 1.28   0.95   6.15 1.01      773      844
##    y1[10]   2.56   2.15 1.66 1.07   0.89   5.50 1.00     1818     1686
##    y1[11]   4.25   3.69 2.50 1.84   1.54   8.59 1.00     1788     1845
##    y1[12]   1.71   1.45 1.09 0.77   0.56   3.62 1.00     1409     1283
##    y1[13]   4.12   3.60 2.34 1.79   1.46   8.59 1.00     1750     1605
##    y1[14]   2.04   1.79 1.14 0.93   0.70   4.16 1.00     1909     1891
##    y1[15]   2.05   1.77 1.28 0.90   0.72   4.26 1.00     1785     1654
##    y1[16]   3.51   3.05 2.05 1.51   1.24   7.14 1.00     1920     1739
##    y1[17]   3.72   3.25 2.16 1.62   1.31   7.45 1.00     1927     1930
##    y1[18]   1.76   1.51 1.08 0.76   0.64   3.77 1.00     1892     1880
##    y1[19]   1.46   1.25 0.91 0.65   0.48   3.10 1.00     1573     1992
##    y1[20]   4.38   3.80 2.91 1.92   1.52   8.94 1.00     1999     2001
##    m1[1]    0.86   0.86 0.18 0.17   0.58   1.13 1.00     1075     1204
##    m1[2]    1.34   1.35 0.23 0.21   0.96   1.71 1.01      805      971
##    m1[3]    1.00   1.00 0.19 0.19   0.69   1.30 1.01      859     1132
##    m1[4]    0.71   0.71 0.20 0.18   0.39   1.02 1.00     1100     1234
##    m1[5]    0.89   0.89 0.17 0.16   0.62   1.15 1.00     1189     1292
##    m1[6]    0.72   0.72 0.18 0.18   0.42   1.01 1.02      336      641
##    m1[7]    0.99   1.00 0.13 0.13   0.77   1.20 1.00     2079     1304
##    m1[8]    1.14   1.15 0.16 0.15   0.87   1.40 1.00     1679     1030
##    m1[9]    0.88   0.87 0.24 0.24   0.48   1.27 1.02      288      298
##    m1[10]   0.77   0.76 0.21 0.21   0.42   1.10 1.01      809     1127
##    m1[11]   1.31   1.31 0.21 0.20   0.95   1.65 1.01      848      888
##    m1[12]   0.37   0.38 0.25 0.24  -0.05   0.76 1.01      516      676
##    m1[13]   1.28   1.28 0.18 0.17   0.99   1.56 1.00     1694     1417
##    m1[14]   0.57   0.58 0.17 0.17   0.29   0.85 1.01      580      987
##    m1[15]   0.59   0.59 0.17 0.16   0.31   0.86 1.01      584      935
##    m1[16]   1.10   1.10 0.18 0.17   0.82   1.39 1.00     1176     1155
##    m1[17]   1.18   1.18 0.17 0.16   0.90   1.45 1.00      929     1034
##    m1[18]   0.41   0.41 0.21 0.19   0.07   0.73 1.00     1844     1032
##    m1[19]   0.21   0.22 0.27 0.24  -0.24   0.62 1.00     1187      966
##    m1[20]   1.32   1.33 0.21 0.19   0.98   1.66 1.00     1036      941

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

poisson regression

objective variable [0,Infinity), integer. variance of error is near to mean  
(normal linear regression, correlation is near to 1,-1, variance of error is constant)  

# of samples is N,  
log li=sum(bj*xji),j=[0,K],i=[1,N]  
yi~Po(li)  
 or  
li=sum(bj*xji),j=[0,k]  
yi~Po(exp li)  

when xj -> xj +1, y -> y* exp bj   

ex6-1.stan

poisson regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] l=X*b;
  y~poisson_log(l);
}
generated quantities{
  array[N] int y1;
  vector[N] l1=X*b;
  for(i in 1:N){
    y1[i]=poisson_log_rng(l1[i]);
  }
}
n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
          y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='poisson')
## 
## Call:  glm(formula = f, family = "poisson", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##       0.314        1.136        0.122  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       57.1 
## Residual Deviance: 37.2  AIC: 104
X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -23.2474 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12       -9.7876   9.45541e-05   0.000384669           1           1       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -9.79
##    b[1]       0.31
##    b[2]       1.14
##    b[3]       0.12
##    y1[1]      1.00
##    y1[2]      3.00
##    y1[3]      1.00
##    y1[4]      2.00
##    y1[5]     12.00
##    y1[6]      0.00
##    y1[7]      3.00
##    y1[8]      1.00
##    y1[9]      3.00
##    y1[10]     2.00
##    y1[11]     1.00
##    y1[12]     3.00
##    y1[13]     0.00
##    y1[14]     2.00
##    y1[15]     3.00
##    y1[16]     0.00
##    y1[17]     0.00
##    y1[18]     6.00
##    y1[19]     0.00
##    y1[20]     2.00
##    y1[21]     1.00
##    y1[22]     1.00
##    y1[23]     3.00
##    y1[24]     2.00
##    y1[25]     1.00
##    y1[26]     3.00
##    y1[27]     4.00
##    y1[28]     1.00
##    y1[29]     2.00
##    y1[30]     1.00
##    l1[1]      0.81
##    l1[2]      0.87
##    l1[3]      0.75
##    l1[4]      0.60
##    l1[5]      1.41
##    l1[6]      1.04
##    l1[7]      0.75
##    l1[8]      0.30
##    l1[9]      0.38
##    l1[10]    -0.54
##    l1[11]     1.20
##    l1[12]     1.45
##    l1[13]     0.38
##    l1[14]    -0.35
##    l1[15]    -0.05
##    l1[16]    -0.65
##    l1[17]     0.55
##    l1[18]     1.42
##    l1[19]    -0.17
##    l1[20]     0.58
##    l1[21]    -0.13
##    l1[22]    -0.36
##    l1[23]     0.81
##    l1[24]     1.00
##    l1[25]     0.00
##    l1[26]     0.93
##    l1[27]     1.44
##    l1[28]     0.48
##    l1[29]    -0.13
##    l1[30]    -0.69
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='l1',ch=F)
##  variable   mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -11.28 -10.98 1.19 0.99 -13.62 -9.97 1.00      855     1191
##    b[1]     0.28   0.28 0.23 0.23  -0.12  0.65 1.00     1095     1155
##    b[2]     1.14   1.13 0.27 0.26   0.71  1.58 1.00     1183     1108
##    b[3]     0.14   0.13 0.27 0.27  -0.30  0.58 1.00      988     1021
##    y1[1]    2.22   2.00 1.51 1.48   0.00  5.00 1.00     1810     1893
##    y1[2]    2.41   2.00 1.64 1.48   0.00  5.00 1.00     1940     1932
##    y1[3]    2.07   2.00 1.47 1.48   0.00  5.00 1.00     2094     1800
##    y1[4]    1.79   2.00 1.36 1.48   0.00  4.00 1.00     1902     1875
##    y1[5]    4.15   4.00 2.19 2.97   1.00  8.00 1.00     2026     1944
##    y1[6]    2.81   3.00 1.76 1.48   0.00  6.00 1.00     2125     2023
##    y1[7]    2.12   2.00 1.55 1.48   0.00  5.00 1.00     1724     1854
##    y1[8]    1.31   1.00 1.21 1.48   0.00  4.00 1.00     1872     1959
##    y1[9]    1.41   1.00 1.21 1.48   0.00  4.00 1.00     2076     2045
##    y1[10]   0.57   0.00 0.79 0.00   0.00  2.00 1.00     2004     1884
##    y1[11]   3.29   3.00 2.03 1.48   0.00  7.00 1.00     1853     1677
##    y1[12]   4.26   4.00 2.26 1.48   1.00  8.00 1.00     2048     2016
##    y1[13]   1.50   1.00 1.25 1.48   0.00  4.00 1.00     2018     2097
##    y1[14]   0.74   1.00 0.88 1.48   0.00  2.00 1.00     1950     1999
##    y1[15]   0.96   1.00 1.03 1.48   0.00  3.00 1.00     1357     1745
##    y1[16]   0.53   0.00 0.77 0.00   0.00  2.00 1.00     1982     1767
##    y1[17]   1.70   2.00 1.35 1.48   0.00  4.00 1.00     1951     1777
##    y1[18]   4.22   4.00 2.35 2.97   1.00  8.00 1.00     1889     1894
##    y1[19]   0.89   1.00 0.96 1.48   0.00  3.00 1.00     1881     1913
##    y1[20]   1.76   2.00 1.34 1.48   0.00  4.00 1.00     2044     2029
##    y1[21]   0.87   1.00 0.95 1.48   0.00  3.00 1.00     1791     1921
##    y1[22]   0.70   0.00 0.86 0.00   0.00  2.00 1.00     2008     1920
##    y1[23]   2.17   2.00 1.51 1.48   0.00  5.00 1.00     2120     1946
##    y1[24]   2.74   3.00 1.66 1.48   0.00  6.00 1.00     1951     1959
##    y1[25]   0.99   1.00 1.05 1.48   0.00  3.00 1.00     1891     2086
##    y1[26]   2.52   2.00 1.67 1.48   0.00  6.00 1.00     1988     1833
##    y1[27]   4.22   4.00 2.29 1.48   1.00  8.00 1.00     1848     1702
##    y1[28]   1.65   1.00 1.34 1.48   0.00  4.00 1.00     1837     1976
##    y1[29]   0.90   1.00 0.97 1.48   0.00  3.00 1.00     2218     1776
##    y1[30]   0.49   0.00 0.76 0.00   0.00  2.00 1.00     1892     1869
##    l1[1]    0.79   0.80 0.17 0.17   0.50  1.07 1.00     2073     1541
##    l1[2]    0.83   0.84 0.22 0.22   0.46  1.17 1.00     1280     1378
##    l1[3]    0.73   0.74 0.17 0.17   0.44  1.01 1.00     2055     1486
##    l1[4]    0.58   0.58 0.18 0.18   0.27  0.87 1.00     1941     1584
##    l1[5]    1.39   1.39 0.21 0.21   1.03  1.73 1.00     1620     1442
##    l1[6]    1.02   1.02 0.18 0.17   0.72  1.30 1.00     1999     1432
##    l1[7]    0.71   0.72 0.21 0.22   0.35  1.04 1.00     1236     1359
##    l1[8]    0.26   0.27 0.24 0.23  -0.15  0.63 1.00     1091     1128
##    l1[9]    0.35   0.36 0.21 0.21   0.00  0.69 1.00     1729     1658
##    l1[10]  -0.57  -0.55 0.38 0.37  -1.20  0.03 1.00     1300     1177
##    l1[11]   1.17   1.18 0.24 0.24   0.76  1.55 1.00     1382     1287
##    l1[12]   1.43   1.43 0.22 0.22   1.07  1.78 1.00     1587     1380
##    l1[13]   0.36   0.36 0.21 0.21   0.01  0.69 1.00     1733     1658
##    l1[14]  -0.39  -0.38 0.33 0.33  -0.95  0.15 1.00     1053     1035
##    l1[15]  -0.09  -0.08 0.28 0.27  -0.57  0.36 1.00     1056     1030
##    l1[16]  -0.69  -0.67 0.39 0.39  -1.33 -0.06 1.00     1057     1038
##    l1[17]   0.52   0.53 0.19 0.19   0.21  0.83 1.00     1895     1659
##    l1[18]   1.38   1.40 0.27 0.27   0.94  1.82 1.00     1412     1259
##    l1[19]  -0.19  -0.18 0.30 0.30  -0.69  0.29 1.00     1387     1311
##    l1[20]   0.55   0.56 0.19 0.19   0.24  0.85 1.00     1921     1611
##    l1[21]  -0.18  -0.16 0.29 0.28  -0.68  0.30 1.00     1052     1042
##    l1[22]  -0.39  -0.37 0.34 0.33  -0.96  0.15 1.00     1330     1258
##    l1[23]   0.78   0.79 0.17 0.17   0.50  1.06 1.00     2072     1542
##    l1[24]   0.98   0.98 0.17 0.17   0.68  1.26 1.00     2031     1433
##    l1[25]  -0.04  -0.03 0.27 0.26  -0.51  0.40 1.00     1058     1041
##    l1[26]   0.90   0.91 0.22 0.22   0.53  1.24 1.00     1308     1377
##    l1[27]   1.42   1.43 0.22 0.22   1.06  1.77 1.00     1592     1380
##    l1[28]   0.45   0.45 0.22 0.22   0.07  0.78 1.00     1139     1248
##    l1[29]  -0.15  -0.14 0.29 0.29  -0.64  0.32 1.00     1402     1409
##    l1[30]  -0.73  -0.71 0.39 0.40  -1.38 -0.09 1.00     1057     1051

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 2 3 4
##   0 3 2 2 0 0
##   1 0 6 2 0 1
##   2 1 1 3 2 0
##   4 0 0 2 1 0
##   5 0 1 0 0 1
##   6 0 0 0 0 1
##   7 0 0 0 0 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 2 3 4
##   0 2 0 1 0 0
##   1 0 5 0 0 0
##   2 0 1 1 0 0
##   4 0 0 1 1 0
##   5 0 0 0 0 1
##   6 0 0 0 0 0
##   7 0 0 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 2 3 4
##   0 1 2 1 0 0
##   1 0 1 2 0 1
##   2 1 0 2 2 0
##   4 0 0 1 0 0
##   5 0 1 0 0 0
##   6 0 0 0 0 1
##   7 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3 4
##   0 4 2 1 0 0
##   1 4 3 1 0 1
##   2 1 3 2 1 0
##   4 0 0 2 1 0
##   5 1 0 0 0 1
##   6 0 0 0 0 1
##   7 0 0 0 1 0
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3 4
##   0 2 0 1 0 0
##   1 4 1 0 0 0
##   2 0 1 1 0 0
##   4 0 0 1 1 0
##   5 0 0 0 0 1
##   6 0 0 0 0 0
##   7 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3 4
##   0 2 2 0 0 0
##   1 0 2 1 0 1
##   2 1 2 1 1 0
##   4 0 0 1 0 0
##   5 1 0 0 0 0
##   6 0 0 0 0 1
##   7 0 0 0 1 0
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

logistic regression

# of samples is N,  
objective variable 0/1 binary  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~Ber(pi), 0/1 binary  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-2.stan

logistic regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~bernoulli_logit(z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=bernoulli_rng(inv_logit(z1[i]));
  }
}
n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='binomial') # it can caluculte when all trials are once
## 
## Call:  glm(formula = f, family = "binomial", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##        0.33         1.03         1.46  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       32.6 
## Residual Deviance: 30.1  AIC: 36.1
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -33.4797 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -15.0554   0.000209731   0.000287005      0.7124      0.7124       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -15.06
##    b[1]       0.33
##    b[2]       1.03
##    b[3]       1.46
##    y1[1]      0.00
##    y1[2]      1.00
##    y1[3]      0.00
##    y1[4]      1.00
##    y1[5]      1.00
##    y1[6]      0.00
##    y1[7]      1.00
##    y1[8]      0.00
##    y1[9]      1.00
##    y1[10]     0.00
##    y1[11]     1.00
##    y1[12]     0.00
##    y1[13]     1.00
##    y1[14]     0.00
##    y1[15]     1.00
##    y1[16]     1.00
##    y1[17]     1.00
##    y1[18]     0.00
##    y1[19]     1.00
##    y1[20]     0.00
##    y1[21]     1.00
##    y1[22]     0.00
##    y1[23]     1.00
##    y1[24]     1.00
##    y1[25]     1.00
##    y1[26]     1.00
##    y1[27]     1.00
##    y1[28]     0.00
##    y1[29]     1.00
##    y1[30]     0.00
##    z1[1]      1.28
##    z1[2]      0.57
##    z1[3]      1.16
##    z1[4]      2.55
##    z1[5]      2.33
##    z1[6]      1.72
##    z1[7]      0.91
##    z1[8]      0.51
##    z1[9]      0.68
##    z1[10]     1.29
##    z1[11]     1.39
##    z1[12]     0.63
##    z1[13]     2.36
##    z1[14]     1.08
##    z1[15]     2.07
##    z1[16]     0.80
##    z1[17]     1.17
##    z1[18]     1.18
##    z1[19]     2.67
##    z1[20]     2.76
##    z1[21]     1.24
##    z1[22]     1.27
##    z1[23]     2.20
##    z1[24]     1.08
##    z1[25]     1.01
##    z1[26]     1.22
##    z1[27]     0.52
##    z1[28]     0.68
##    z1[29]     1.97
##    z1[30]    -0.54
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.74 -16.36 1.39 1.11 -19.50 -15.24 1.00      678      776
##    b[1]     0.28   0.26 0.77 0.73  -0.96   1.55 1.00      873      916
##    b[2]     1.29   1.24 1.11 1.07  -0.41   3.23 1.01      786      551
##    b[3]     1.86   1.78 1.27 1.17  -0.13   4.04 1.00      697      628
##    y1[1]    0.79   1.00 0.41 0.00   0.00   1.00 1.00     1653       NA
##    y1[2]    0.63   1.00 0.48 0.00   0.00   1.00 1.00     1912       NA
##    y1[3]    0.77   1.00 0.42 0.00   0.00   1.00 1.00     2112       NA
##    y1[4]    0.92   1.00 0.27 0.00   0.00   1.00 1.00     1780       NA
##    y1[5]    0.91   1.00 0.28 0.00   0.00   1.00 1.00     1790       NA
##    y1[6]    0.84   1.00 0.36 0.00   0.00   1.00 1.00     1941       NA
##    y1[7]    0.74   1.00 0.44 0.00   0.00   1.00 1.00     2054       NA
##    y1[8]    0.60   1.00 0.49 0.00   0.00   1.00 1.00     1793       NA
##    y1[9]    0.65   1.00 0.48 0.00   0.00   1.00 1.00     1996       NA
##    y1[10]   0.79   1.00 0.41 0.00   0.00   1.00 1.00     2002       NA
##    y1[11]   0.81   1.00 0.40 0.00   0.00   1.00 1.00     1942       NA
##    y1[12]   0.66   1.00 0.47 0.00   0.00   1.00 1.00     1799       NA
##    y1[13]   0.92   1.00 0.28 0.00   0.00   1.00 1.00     1630       NA
##    y1[14]   0.75   1.00 0.43 0.00   0.00   1.00 1.00     1824       NA
##    y1[15]   0.90   1.00 0.30 0.00   0.00   1.00 1.00     1837       NA
##    y1[16]   0.68   1.00 0.47 0.00   0.00   1.00 1.00     1926       NA
##    y1[17]   0.76   1.00 0.42 0.00   0.00   1.00 1.00     2035       NA
##    y1[18]   0.78   1.00 0.42 0.00   0.00   1.00 1.00     1997       NA
##    y1[19]   0.92   1.00 0.27 0.00   0.00   1.00 1.00     1999       NA
##    y1[20]   0.93   1.00 0.26 0.00   0.00   1.00 1.00     1484       NA
##    y1[21]   0.79   1.00 0.41 0.00   0.00   1.00 1.00     1874       NA
##    y1[22]   0.80   1.00 0.40 0.00   0.00   1.00 1.00     1998       NA
##    y1[23]   0.91   1.00 0.28 0.00   0.00   1.00 1.00     1790       NA
##    y1[24]   0.76   1.00 0.43 0.00   0.00   1.00 1.00     1738       NA
##    y1[25]   0.74   1.00 0.44 0.00   0.00   1.00 1.00     1789       NA
##    y1[26]   0.76   1.00 0.43 0.00   0.00   1.00 1.00     2004       NA
##    y1[27]   0.60   1.00 0.49 0.00   0.00   1.00 1.00     1878       NA
##    y1[28]   0.67   1.00 0.47 0.00   0.00   1.00 1.00     1897       NA
##    y1[29]   0.89   1.00 0.32 0.00   0.00   1.00 1.00     1697       NA
##    y1[30]   0.36   0.00 0.48 0.00   0.00   1.00 1.00     1518       NA
##    z1[1]    1.47   1.42 0.80 0.79   0.24   2.86 1.00     1555     1475
##    z1[2]    0.58   0.57 0.64 0.63  -0.43   1.65 1.00     1171     1161
##    z1[3]    1.35   1.28 1.06 0.96  -0.32   3.17 1.00     1671     1286
##    z1[4]    3.08   2.97 1.33 1.29   1.17   5.49 1.00      696      728
##    z1[5]    2.81   2.70 1.17 1.12   1.11   4.90 1.00      742      746
##    z1[6]    2.05   1.98 0.92 0.83   0.69   3.66 1.00     1216     1259
##    z1[7]    1.01   0.97 0.61 0.60   0.06   2.08 1.00     1830     1483
##    z1[8]    0.50   0.49 0.67 0.65  -0.55   1.63 1.00     1070     1044
##    z1[9]    0.71   0.70 0.61 0.60  -0.23   1.73 1.00     1358     1296
##    z1[10]   1.48   1.43 0.80 0.80   0.25   2.87 1.00     1544     1474
##    z1[11]   1.64   1.56 0.96 0.85   0.20   3.34 1.00     1692     1364
##    z1[12]   0.66   0.65 0.62 0.60  -0.31   1.70 1.00     1274     1284
##    z1[13]   2.85   2.74 1.19 1.15   1.13   4.97 1.00      733      761
##    z1[14]   1.22   1.18 0.68 0.66   0.17   2.42 1.00     1796     1373
##    z1[15]   2.49   2.39 1.02 0.95   1.03   4.34 1.00      842      806
##    z1[16]   0.86   0.84 0.60 0.59  -0.07   1.87 1.00     1639     1345
##    z1[17]   1.33   1.30 0.73 0.72   0.22   2.60 1.00     1688     1483
##    z1[18]   1.34   1.31 0.73 0.73   0.22   2.62 1.00     1674     1506
##    z1[19]   3.24   3.12 1.43 1.37   1.15   5.80 1.00      683      722
##    z1[20]   3.35   3.23 1.51 1.44   1.14   6.05 1.00      674      705
##    z1[21]   1.45   1.37 1.01 0.93  -0.13   3.23 1.00     1695     1341
##    z1[22]   1.46   1.42 0.79 0.79   0.24   2.84 1.00     1562     1475
##    z1[23]   2.65   2.55 1.09 1.03   1.11   4.62 1.00      785      802
##    z1[24]   1.24   1.17 1.10 1.01  -0.50   3.14 1.00     1639     1220
##    z1[25]   1.13   1.10 0.64 0.62   0.13   2.26 1.00     1847     1485
##    z1[26]   1.42   1.34 1.03 0.93  -0.19   3.20 1.00     1689     1304
##    z1[27]   0.52   0.51 0.66 0.65  -0.52   1.63 1.00     1098     1090
##    z1[28]   0.72   0.71 0.60 0.59  -0.22   1.74 1.00     1374     1314
##    z1[29]   2.37   2.26 0.98 0.90   0.96   4.08 1.00      912      941
##    z1[30]  -0.82  -0.75 1.55 1.48  -3.42   1.64 1.01      718      639

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##      0  1
##   0  1  6
##   1  0 23
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##      0  1
##   0  1  4
##   1  0 11
## 
## , ,  = b
## 
##    
##      0  1
##   0  0  2
##   1  0 12
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##      0  1
##   0  1  6
##   1  0 23
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##      0  1
##   0  1  4
##   1  0 11
## 
## , ,  = b
## 
##    map
##      0  1
##   0  0  2
##   1  0 12
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

multi logistic regression

# of samples is N,  
# of trials of a sample i is mi,  
objective variable [0,n], integer  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~B(mi, pi), # of acutual incident  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-3.stan

multi logistic regression

data{
  int N;
  int K;
  array[N] int m;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~binomial_logit(m,z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
  }
}
n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)

mdl=cmdstan_model('./ex6-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -131.051 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -87.7291    0.00014274   0.000204344      0.8569      0.8569       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -87.73
##    b[1]       0.34
##    b[2]       1.39
##    b[3]       1.01
##    y1[1]      5.00
##    y1[2]      8.00
##    y1[3]      1.00
##    y1[4]      7.00
##    y1[5]      0.00
##    y1[6]      3.00
##    y1[7]      7.00
##    y1[8]      2.00
##    y1[9]      1.00
##    y1[10]     7.00
##    y1[11]     0.00
##    y1[12]     5.00
##    y1[13]     2.00
##    y1[14]     2.00
##    y1[15]     3.00
##    y1[16]     1.00
##    y1[17]     3.00
##    y1[18]     6.00
##    y1[19]     6.00
##    y1[20]     3.00
##    y1[21]     4.00
##    y1[22]     6.00
##    y1[23]     8.00
##    y1[24]     0.00
##    y1[25]     1.00
##    y1[26]     7.00
##    y1[27]     3.00
##    y1[28]     3.00
##    y1[29]     6.00
##    y1[30]     2.00
##    z1[1]      0.47
##    z1[2]      1.08
##    z1[3]     -0.10
##    z1[4]      1.90
##    z1[5]     -0.09
##    z1[6]      0.06
##    z1[7]      2.49
##    z1[8]      0.21
##    z1[9]      0.36
##    z1[10]     0.71
##    z1[11]     0.84
##    z1[12]     2.42
##    z1[13]    -0.70
##    z1[14]     1.93
##    z1[15]     0.90
##    z1[16]     1.05
##    z1[17]    -0.20
##    z1[18]     2.67
##    z1[19]     0.45
##    z1[20]     1.14
##    z1[21]     0.08
##    z1[22]     2.68
##    z1[23]     1.99
##    z1[24]     0.27
##    z1[25]     2.16
##    z1[26]     1.85
##    z1[27]    -0.56
##    z1[28]     0.12
##    z1[29]     1.02
##    z1[30]     0.41
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -89.39 -88.98 1.39 1.06 -92.11 -87.94 1.00      820      864
##    b[1]     0.33   0.33 0.34 0.33  -0.24   0.88 1.00      670      595
##    b[2]     1.44   1.43 0.40 0.39   0.83   2.11 1.00     1115      930
##    b[3]     1.07   1.06 0.41 0.39   0.43   1.78 1.00      720      693
##    y1[1]    4.26   4.00 1.38 1.48   2.00   6.00 1.00     2010     1888
##    y1[2]    6.74   7.00 1.38 1.48   4.00   9.00 1.00     2086       NA
##    y1[3]    3.26   3.00 1.41 1.48   1.00   6.00 1.00     1826     1859
##    y1[4]    6.12   6.00 0.90 1.48   4.00   7.00 1.00     1993       NA
##    y1[5]    0.47   0.00 0.50 0.00   0.00   1.00 1.00     2052       NA
##    y1[6]    3.06   3.00 1.29 1.48   1.00   5.00 1.00     1915     2073
##    y1[7]    6.46   7.00 0.72 0.00   5.00   7.00 1.00     1986       NA
##    y1[8]    3.87   4.00 1.43 1.48   1.00   6.00 1.00     1837     1933
##    y1[9]    0.59   1.00 0.49 0.00   0.00   1.00 1.00     1999       NA
##    y1[10]   5.38   5.00 1.41 1.48   3.00   8.00 1.00     1965       NA
##    y1[11]   0.69   1.00 0.46 0.00   0.00   1.00 1.00     1898       NA
##    y1[12]   5.51   6.00 0.69 0.00   4.00   6.00 1.00     1953       NA
##    y1[13]   2.62   2.00 1.46 1.48   0.00   5.00 1.00     1883     1834
##    y1[14]   1.75   2.00 0.47 0.00   1.00   2.00 1.00     2025       NA
##    y1[15]   4.98   5.00 1.24 1.48   3.00   7.00 1.00     2015       NA
##    y1[16]   0.73   1.00 0.44 0.00   0.00   1.00 1.00     1765       NA
##    y1[17]   2.63   3.00 1.28 1.48   1.00   5.00 1.00     1718     1788
##    y1[18]   5.58   6.00 0.65 0.00   4.00   6.00 1.00     1836       NA
##    y1[19]   5.45   5.00 1.57 1.48   3.00   8.00 1.00     1482     1745
##    y1[20]   2.26   2.00 0.77 1.48   1.00   3.00 1.00     1849       NA
##    y1[21]   2.60   3.00 1.18 1.48   1.00   4.00 1.00     1694     1588
##    y1[22]   7.46   8.00 0.73 0.00   6.00   8.00 1.00     1661       NA
##    y1[23]   7.96   8.00 1.02 1.48   6.00   9.00 1.00     1977       NA
##    y1[24]   0.57   1.00 0.49 0.00   0.00   1.00 1.00     1922       NA
##    y1[25]   0.89   1.00 0.31 0.00   0.00   1.00 1.00     2024       NA
##    y1[26]   6.96   7.00 1.00 1.48   5.00   8.00 1.00     1842       NA
##    y1[27]   1.43   1.00 1.00 1.48   0.00   3.00 1.00     1859     1781
##    y1[28]   3.18   3.00 1.31 1.48   1.00   5.00 1.00     1900     1929
##    y1[29]   5.93   6.00 1.28 1.48   4.00   8.00 1.00     1935       NA
##    y1[30]   2.38   2.00 0.98 1.48   1.00   4.00 1.00     1965       NA
##    z1[1]    0.48   0.47 0.28 0.27   0.03   0.95 1.00     1451     1682
##    z1[2]    1.11   1.10 0.24 0.25   0.74   1.51 1.00     1716     1562
##    z1[3]   -0.14  -0.12 0.33 0.32  -0.70   0.37 1.00      782      731
##    z1[4]    1.97   1.95 0.36 0.36   1.40   2.57 1.00     1444     1276
##    z1[5]   -0.12  -0.11 0.33 0.32  -0.69   0.38 1.00      777      724
##    z1[6]    0.05   0.04 0.35 0.36  -0.52   0.64 1.00     1277     1598
##    z1[7]    2.58   2.56 0.50 0.48   1.81   3.42 1.01     1303     1174
##    z1[8]    0.21   0.20 0.32 0.33  -0.31   0.74 1.00     1332     1637
##    z1[9]    0.34   0.35 0.34 0.33  -0.22   0.90 1.00      669      595
##    z1[10]   0.73   0.73 0.25 0.25   0.33   1.16 1.00     1578     1670
##    z1[11]   0.84   0.84 0.41 0.41   0.16   1.52 1.00      706      609
##    z1[12]   2.51   2.49 0.48 0.47   1.77   3.32 1.01     1316     1189
##    z1[13]  -0.76  -0.75 0.38 0.37  -1.41  -0.16 1.00     1142      904
##    z1[14]   1.99   1.97 0.36 0.36   1.42   2.60 1.00     1436     1243
##    z1[15]   0.93   0.92 0.24 0.24   0.55   1.33 1.00     1672     1614
##    z1[16]   1.06   1.05 0.45 0.44   0.32   1.81 1.00      756      833
##    z1[17]  -0.24  -0.23 0.33 0.32  -0.81   0.28 1.00      831      815
##    z1[18]   2.77   2.74 0.54 0.53   1.93   3.69 1.01     1276     1158
##    z1[19]   0.44   0.44 0.35 0.34  -0.15   1.00 1.00      662      600
##    z1[20]   1.15   1.14 0.46 0.46   0.40   1.93 1.00      781      807
##    z1[21]   0.05   0.06 0.33 0.32  -0.51   0.57 1.00      720      579
##    z1[22]   2.78   2.75 0.54 0.53   1.94   3.70 1.01     1275     1158
##    z1[23]   2.06   2.04 0.38 0.38   1.47   2.69 1.00     1415     1217
##    z1[24]   0.27   0.26 0.31 0.31  -0.24   0.79 1.00     1355     1652
##    z1[25]   2.24   2.22 0.42 0.41   1.59   2.95 1.01     1369     1165
##    z1[26]   1.92   1.90 0.35 0.35   1.37   2.50 1.00     1457     1374
##    z1[27]  -0.61  -0.59 0.36 0.35  -1.22  -0.05 1.00     1100      863
##    z1[28]   0.12   0.11 0.34 0.34  -0.43   0.68 1.00     1297     1571
##    z1[29]   1.05   1.05 0.24 0.24   0.68   1.45 1.00     1702     1502
##    z1[30]   0.41   0.41 0.29 0.28  -0.05   0.90 1.00     1427     1713

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 2 3 4 5 6 7 8
##   0 1 0 0 0 0 0 0 0 0
##   1 0 6 0 0 0 0 0 0 0
##   2 0 0 1 3 0 0 0 0 0
##   3 0 0 2 1 1 0 0 0 0
##   4 0 0 1 1 0 1 0 0 0
##   5 0 0 0 0 1 1 1 0 0
##   6 0 0 0 0 0 1 3 1 0
##   7 0 0 0 0 0 0 0 1 2
##   8 0 0 0 0 0 0 0 1 0
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 2 3 4 5 6 7 8
##   0 1 0 0 0 0 0 0 0 0
##   1 0 4 0 0 0 0 0 0 0
##   2 0 0 0 2 0 0 0 0 0
##   3 0 0 1 1 0 0 0 0 0
##   4 0 0 1 0 0 0 0 0 0
##   5 0 0 0 0 0 1 0 0 0
##   6 0 0 0 0 0 0 0 0 0
##   7 0 0 0 0 0 0 0 0 0
##   8 0 0 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 2 3 4 5 6 7 8
##   0 0 0 0 0 0 0 0 0 0
##   1 0 2 0 0 0 0 0 0 0
##   2 0 0 1 1 0 0 0 0 0
##   3 0 0 1 0 1 0 0 0 0
##   4 0 0 0 1 0 1 0 0 0
##   5 0 0 0 0 1 0 1 0 0
##   6 0 0 0 0 0 1 3 1 0
##   7 0 0 0 0 0 0 0 1 2
##   8 0 0 0 0 0 0 0 1 0
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3 4 5 6 7 8
##   0 1 0 0 0 0 0 0 0 0
##   1 0 6 0 0 0 0 0 0 0
##   2 0 0 1 3 0 0 0 0 0
##   3 0 0 1 2 1 0 0 0 0
##   4 0 0 1 1 0 1 0 0 0
##   5 0 0 0 0 1 1 1 0 0
##   6 0 0 0 0 0 0 4 1 0
##   7 0 0 0 0 0 0 0 1 2
##   8 0 0 0 0 0 0 0 1 0
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3 4 5 6 7 8
##   0 1 0 0 0 0 0 0 0 0
##   1 0 4 0 0 0 0 0 0 0
##   2 0 0 0 2 0 0 0 0 0
##   3 0 0 0 2 0 0 0 0 0
##   4 0 0 1 0 0 0 0 0 0
##   5 0 0 0 0 0 1 0 0 0
##   6 0 0 0 0 0 0 0 0 0
##   7 0 0 0 0 0 0 0 0 0
##   8 0 0 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3 4 5 6 7 8
##   0 0 0 0 0 0 0 0 0 0
##   1 0 2 0 0 0 0 0 0 0
##   2 0 0 1 1 0 0 0 0 0
##   3 0 0 1 0 1 0 0 0 0
##   4 0 0 0 1 0 1 0 0 0
##   5 0 0 0 0 1 0 1 0 0
##   6 0 0 0 0 0 0 4 1 0
##   7 0 0 0 0 0 0 0 1 2
##   8 0 0 0 0 0 0 0 1 0
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

gamma regression

objective variable [0,Infinity)

# of samples is N,  
log (a/li)=sum(bj*xji),j=[0,K],i=[1,N]
li=a/exp(sum(bj*xji))
yi~Ga(a,li)  

ex6-4.stan

gamma regression

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> a;
}
model{
  vector[N] l;
  for(i in 1:N){
    l[i]=a/exp(X[i]*b);
  }
  y~gamma(a,l);
}
n=20
tb=tibble(x1=runif(n,0,2),x2=runif(n,0,2),
          y=rgamma(n,3,3/exp(x1+x2)))

f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-4.stan')

mle=mdl$optimize(data=data)  # it sometimes occur error and stop process
## Initial log joint probability = -1039.58 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -52.3623   5.79493e-05    0.00064872      0.5011           1       21    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -52.36
##      b[1]    -1.03
##      b[2]     1.59
##      b[3]     1.41
##      a        3.95
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -53.12 -52.74 1.60 1.26 -56.31 -51.33 1.00      618      530
##      b[1]  -1.00  -1.00 0.40 0.37  -1.63  -0.35 1.00      656      596
##      b[2]   1.59   1.58 0.21 0.20   1.24   1.94 1.00      855      902
##      b[3]   1.40   1.40 0.25 0.24   1.02   1.80 1.00      845      856
##      a      3.82   3.66 1.14 1.06   2.23   5.99 1.00      945     1217

negative binomial regression

The event with probability p do not occur y times till the event occur a times  
(negative binomial distribution has larger variance than poisson distribution)

y~NB(a,p), log(p)=X*b

y~NB(a,l0/(1+l0)) when y~Po(l), l~Ga(a,l0), l0=a/E[l]

ex6-5.stan

negative binomial regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> a;
}
model{
  a~cauchy(0,5);
  y~neg_binomial_2_log(X*b,a);
}
n=20
tb=tibble(x1=runif(n,-1,0),x2=runif(n,-1,0),
          y=rnbinom(n,3,exp(x1+x2)))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-5.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -169.256 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -46.4117   0.000303863   6.78782e-05           1           1       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -46.41
##      b[1]    -0.57
##      b[2]    -2.53
##      b[3]    -1.22
##      a        3.89
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -47.05 -46.68 1.56 1.29 -50.05 -45.27 1.00      792     1024
##      b[1]  -0.63  -0.62 0.51 0.49  -1.44   0.17 1.00      903     1054
##      b[2]  -2.60  -2.59 0.60 0.59  -3.57  -1.65 1.00     1029     1045
##      b[3]  -1.27  -1.27 0.55 0.52  -2.19  -0.35 1.01     1053     1023
##      a      5.42   4.45 3.84 2.44   1.74  11.87 1.00     1216      859

beta regression

using prior of binomial distribution parameter p[0,1]
y~B(n,p), p~beta(a,b)

m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)

m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*s

ex6-6.stan

beta regression

data{
  int N;
  int K;
  vector[N] p;
  matrix[N,K] X;
}
parameters{
  vector[K] be;
  real<lower=0> s;
}
model{
  vector[N] a;
  vector[N] b;
  for(i in 1:N){
    a[i]=inv_logit(X[i]*be)*s;
    b[i]=(1-inv_logit(X[i]*be))*s;
  }
  p~beta(a,b);
}
n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
          p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$p)
plot(tb$x2,tb$p)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),p=tb$p,X=X)

mdl=cmdstan_model('./ex6-6.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -66.1095 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23       14.9812    0.00203034   0.000370719           1           1       25    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__     14.98
##     be[1]    -3.28
##     be[2]     6.12
##     be[3]     6.63
##     s         4.85
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##     lp__  14.54  14.86 1.43 1.30 11.73 16.25 1.00      633     1147
##     be[1] -3.35  -3.33 0.60 0.60 -4.36 -2.35 1.00      989      982
##     be[2]  6.27   6.20 1.58 1.57  3.76  8.93 1.00      699      545
##     be[3]  6.73   6.77 1.70 1.60  3.86  9.51 1.01      837      822
##     s      4.95   4.80 1.51 1.46  2.82  7.70 1.00     1418     1481

beta binomial regression

fitting to distribution has larger variance than binomial distribution
y~betaB(n,a,b) when y~B(n,p), p~beta(a,b)

m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)

m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*s

ex6-7.stan

beta binomial regression

data{
  int N;
  int K;
  array[N] int y;
  array[N] int n;
  matrix[N,K] X;
}
parameters{
  vector[K] be;
  real<lower=0> s;
}
model{
  vector[N] a;
  vector[N] b;
  for(i in 1:N){
    a[i]=inv_logit(X[i]*be)*s;
    b[i]=(1-inv_logit(X[i]*be))*s;
  }
  y~beta_binomial(n,a,b);
  s~cauchy(0,5);
}
n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
          p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3),
          n1=floor(runif(n,5,9)),
          y=rbinom(n,n1,p))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),n=tb$n1,y=tb$y,X=X)

mdl=cmdstan_model('./ex6-7.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -89.2412 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -60.8222    0.00598732   0.000460707           1           1       21    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__    -60.82
##     be[1]    -3.83
##     be[2]     3.94
##     be[3]     9.21
##     s         6.08
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##     lp__  -60.93 -60.54   1.54 1.32 -63.98 -59.15 1.01      546      781
##     be[1]  -4.08  -4.08   0.78 0.76  -5.36  -2.79 1.00      856     1006
##     be[2]   4.31   4.31   1.65 1.62   1.52   6.95 1.00      932      611
##     be[3]   9.81   9.73   2.24 2.15   6.12  13.62 1.00      771      799
##     s      49.26  10.21 455.82 7.43   3.39  84.35 1.00      915      545